\documentclass[a4paper, 11pt]{article}
%\usepackage{qtree} %z jakiegoś dziwacznego powodu ten pakiet musi być aż tak wysoko...
\usepackage[OT4,plmath,MeX]{polski}
\usepackage{amstext, amsopn, amsfonts, amsthm, amsbsy, latexsym, amssymb}
\usepackage[utf8]{inputenc}

\usepackage{indentfirst}
\usepackage[margin=2cm]{geometry}
\usepackage{amsmath} %found in http://en.wikibooks.org/wiki/LaTeX/Formatting
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{enumerate}
\usepackage{textcomp}

\usepackage{multirow}

\usepackage{graphicx}
\usepackage{epstopdf}
\usepackage{url}
\usepackage{hyperref}
\usepackage{multirow}
\usepackage{listings}

\newcommand{\HRule}{\rule{\linewidth}{0.5mm}}

\begin{document}

\begin{titlepage}
\begin{center}

% Upper part of the page. The '~' is needed because \\
% only works if a paragraph has started.
\includegraphics[width=0.15\textwidth]{./include/agh-logo}~\\[1cm]

\textsc{\LARGE Akademia Górniczo-Hutnicza w Krakowie}\\
\textsc{\Large Wydział Elektrotechniki, Automatyki, Informatyki i Inżynierii Biomedycznej}\\
[1.5cm]

\textsc{\Large Projekt zaliczeniowy}\\[0.5cm]

% Title
\HRule \\[0.4cm]
{ \huge \bfseries Analizator EKG}\\[0.4cm]

\HRule \\[1.5cm]

% Author and supervisor
\begin{minipage}{0.27\textwidth}
\begin{flushleft} \large
\emph{Autorzy:}\\
 Mateusz \textsc{Baran} \\
 Krzysztof \textsc{Bębenek} \\
 Bartłomiej \textsc{Bułat} \\
 Szczepan \textsc{Czaicki} \\
 Tomasz \textsc{Drzewiecki} \\
 Krzysztof \textsc{Farganus} \\
 Łukasz \textsc{Jaromi} \\
 Mateusz \textsc{Krasucki} \\
 Łukasz \textsc{Krzyżek} \\
 Łukasz \textsc{Kutrzuba} \\
 Weronika \textsc{Łabaj} \\
 Paweł \textsc{Maślanka} \\
\end{flushleft}
\end{minipage}
\begin{minipage}{0.27\textwidth}
\begin{flushleft} \large
 \phantom{Autorzy:}
 Piotr \textsc{Matuszkiewicz} \\
 Norbert \textsc{Pabian} \\
 Łukasz \textsc{Pękala} \\
 Krzysztof \textsc{Piekutowski} \\
 Grzegorz \textsc{Pietrzyk} \\
 Łukasz \textsc{Podolski} \\
 Mikołaj \textsc{Rzepka} \\
 Agata \textsc{Sitnik} \\
 Leszek \textsc{Sosnowski} \\
 Aleksander \textsc{Steliga} \\
 Mateusz \textsc{Ślażyński} \\
 Łukasz \textsc{Zieńkowski}
\end{flushleft}
\end{minipage}
\begin{minipage}{0.4\textwidth}
\begin{flushright} \large
\emph{Opiekun:} \\
 mgr inż. Tomasz \textsc{Pięciak}
\end{flushright}
\end{minipage}

\vfill

% Bottom of the page
{\large \today}

\end{center}
\end{titlepage}


\setcounter{page}{2}
\tableofcontents
\newpage

\section{Specyfikacja zadania}

Celem projektu jest stworzenie zintegrowanego systemu pozwalającego na przeglądanie i automatyczną analizę sygnału EKG. Sygnał dostarczany jest w formie cyfrowej w standardzie wykorzystywanym w MIT–BIH Arrhythmia Database. Różne etapy przetwarzania, takie jak usuwanie linii bazowej, detekcja załamków R czy klasyfikacja zespołów QRS wykonywana jest przez różne zespoły (szczegółowe opisy specyfikacji modułów: \ref{sec:mod}), których praca składa się na jeden program.

Moduły przetwarzania integrowane i uzupełniane są modułami kontrolującymi przepływ danych i odpowiadającymi za komunikację z użytkownikiem. Wzajemne zależności pomiędzy modułami prezentuje rys. \ref{fig:zaleznosci}.

\begin{figure}[h!]
  \centering
  \includegraphics[width=0.7\linewidth]{include/Projekt_zaleznosci}
  \label{fig:zaleznosci}
  \caption{Zależności pomiędzy modułami projektu.}
\end{figure}

\section{Specyfikacja techniczna rozwiązania}
\label{sec:techspec}

\subsection{Wykorzystane narzędzia}
\label{sec:tools}

Podczas realizacji projektu wykorzystywane były różne narzędzia do tworzenia i prototypowania rozwiązań. Wstępne projekty przygotowywane były w programie Matlab, zaś ostateczny kod powstawał w języku C++ (standard '03 z elementami standardu C++11 obsługiwanymi przez wspierane kompilatory).
Minimalne wymagania kompilacji projektu są następujące:
\begin{itemize}
  \item Jeden z kompilatorów:
  \begin{itemize}
    \item Microsoft Visual Studio 2010
    \item GCC 4.5
  \end{itemize}
  \item Biblioteki:
  \begin{itemize}
    \item Boost 1.51
    \item Qt 4.8
    \item Qwt 6.01
    \item gsl 1.15
    \item WFDB 10.5.16 (zawarta w źródłach projektu)
    \item FFTW 3.3.2
   \item KissFFT 1.3.0  (zawarta w źródłach projektu)
   \item ALGLIB 3.7.0   (zawarta w źródłach projektu)
  \end{itemize}
\end{itemize}

Do wersjonowania i śledzenia błędów wykorzystywaliśmy platformę Github wraz z rozproszonym systemem kontroli wersji Git. Posiada on zaawansowane możliwości wspierające pracę grupową nad projektem, co szczególnie przydaje się, gdy liczba osób jest duża.

\subsection{Projekt systemu}
\label{sec:sys_proj}

Program został wykonany w architekturze MVC -- istnieje ścisły podział na część wyświetlającą interfejs użytkownika, moduły przetwarzania sygnału oraz kontroler łączący te dwa elementy, co obrazuje diagram pakietów \ref{fig:package_diagram}. Zastosowano obiektowe podejście przy projektowaniu hierarchii klas realizujących przetwarzanie (rys. \ref{fig:class_diagram}), umożliwiające bezproblemową wymianę implementacji dowolnego modułu na inną.

\begin{figure}[h!]
  \centering
  \includegraphics[width=\linewidth]{include/package_diagram}
  \label{fig:package_diagram}
  \caption{Diagram pakietów.}
\end{figure}

Samo przetwarzanie realizowane jest wielofazowo -- moduły przeliczane są sekwencyjnie na osobnym wątku niż wątek zdarzeń GUI, co pozwala na zatrzymanie zbyt długo trwającej operacji. Zaimplementowano także mechanizm buforowania wyników, dzięki czemu przy zmianie parametrów pewnego modułu nie ma konieczności przeliczania wyników modułów wcześniejszych.

Cała aplikacja napisana jest jako samodzielny, wieloplatformowy program (rys. \ref{fig:deployment_diagram}) -- możliwa jest kompilacja pliku binarnego pod systemami Windows, Linux i MacOS X. Całe przetwarzanie wykonywane jest lokalnie, dzięki czemu nie jest wymagane połączenie z Internetem. Dzięki wybraniu licencji GPL v2+ projekt jest wolny i możliwy jest jego dalszy rozwój i kompilacja na nowe platformy.

\begin{figure}[h!]
  \centering
  \includegraphics[width=0.5\linewidth]{include/deployment_diagram}
  \label{fig:deployment_diagram}
  \caption{Diagram wdrożenia.}
\end{figure}

Diagram sekwencji przedstawiający sposób uruchamiania modułów jest przedstawiony na rys. \ref{fig:sequence_diagram}. Każdy moduł otrzymuje najpierw zestaw parametrów (o ile posiada jakieś parametry), następnie zaś jest uruchamiany funkcją runModule. Kolejne moduły uruchamiane są sekwencyjnie według identycznego schematu. Wskutek akcji użytkownika możliwe jest przerwanie analizy pomiędzy poszczególnymi modułami.

\begin{figure}[h!]
  \centering
  \includegraphics[width=0.5\linewidth]{include/sequence_diagram}
  \caption{Diagram sekwencji.}
  \label{fig:sequence_diagram}
\end{figure}

Diagramy klas zostały przedstawione, z~uwagi na wielkość na trzech rysunkach. Rys. \ref{fig:controller-modules} obrazuje zależności pomiędzy poszczególnymi klasami zrealizowanymi w ramach modułów, klasami abstrakcyjnymi oraz kontrolerem. Rys. \ref{fig:modules-model} przedstawia zależności pomiędzy klasami modułów a~klasami realizującymi model. Rys \ref{fig:controller-gui} przedstawia połączenie kontrolera z klasą odpowiadającą za GUI.

\begin{figure}[h!]
  \centering
  \includegraphics[width=0.9\linewidth]{include/controller-modules}
  \caption{Diagram klas kontrolera i~modułów.}
  \label{fig:controller-modules}
\end{figure}

\begin{figure}[h!]
  \centering
  \includegraphics[width=0.9\linewidth]{include/modules-model}
  \caption{Diagram klas modułów i modelu.}
  \label{fig:modules-model}
\end{figure}

\begin{figure}[h!]
  \centering
  \includegraphics[width=0.9\linewidth]{include/controller-gui}
  \caption{Diagram klas kontrolera i GUI.}
  \label{fig:controller-gui}
\end{figure}

Diagram przypadków użycia, przedstawiony na rys. \ref{fig:use_case_diagram} pokazuje różne scenariusze, w których może być używany zrealizowany system. Każdemu przypadkowi odpowiada jeden lub więcej modułów. zachodzące pomiędzy nimi zależności są ukazywane relacjami zawierania.
\begin{figure}[h!]
  \centering
  \includegraphics[width=\linewidth]{include/use_case_diagram}
  \caption{Diagram klas.}
  \label{fig:use_case_diagram}
\end{figure}

\subsection{Główne elementy kodu systemu}

Głównym elementem systemu jest klasa \verb+ECGController+ \ref{fig:class_diagram}. Metoda \verb+rerunAnalysis+, która się w~niej znajduje odpowiada za kilka rzeczy:
\begin{enumerate}
\item Upewnianie się, że moduły mają aktualne parametry.
\item Rozpoczynanie analizy zapisu EKG.
\item Przerywanie analizy w przypadku, gdy już ona trwa.
\item Raportowanie do GUI o kolejnych etapach pracy systemu (który moduł aktualnie się wykonuje).
\end{enumerate}
Ostatnie zadanie wykonywane jest za pomocą dwóch funkcyjnych argumentów wejściowych.

Parametry, ustawiane metodami \verb+setParam*+ mają postać mapy przypisującej nazwom poszczególnych parametrów odpowiadające im wartości liczbowe. Wykorzystywane są przede wszystkim w~modułach usuwania linii bazowej oraz wykrywania załamków R.

Do klas pomocniczych można zaliczyć między innymi \verb+IntSignal+ czy \verb+ECGSignalChannel+. Sygnały, które przechodzą pomiędzy systemami czynią to właśnie poprzez obiekty tych klas. Są one zrealizowane jako sprytne wskaźniki do struktur trzymających wektory w stylu C biblioteki GSL, dzięki czemu możliwe jest pogodzenie szybkości działania z wygodną składnią i semantyką sprytnego wskaźnika.

Poszczególne moduły są tworzone w hierarchii obiektowej (rys. \ref{fig:class_diagram}). Umożliwia to podmianę w~przyszłości modułu na lepszą implementację bez poprawiania reszty silnika oraz przyczynia się do zmniejszenia współzależności.

\section{Opisy modułów}
\label{sec:mod}

\subsection{Usuwanie linii bazowej}
\label{sec:baseline}
Autorzy: Weronika Łabaj i Piotr Matuszkiewicz.

\subsubsection{Opis zadania}
\label{sec:baseline:desc}

\begin{description}
\item[Temat] Metody filtracji i detekcji izolinii w sygnale EKG
\item[Opis] Występujące w sygnale EKG zakłócenia sieciowe i mięśniowe, jak również falowanie linii izoelektrycznej niejednokrotnie uniemożliwiają przeprowadzenie właściwej i poprawną analizy sygnału. W szczególności mają one wpływ na detekcję odcinków ST. W związku z tym detekcja i usunięcie izolinii, zwanej również linią bazową (\textit{ang. baseline}) jest standardowym pierwszym elementem w~tym procesie. Od tego, w jakim stopniu sygnał oddany do dalszej analizy pozbawiony będzie zakłóceń zależy jej dokładność i powodzenie. Celem niniejszej części projektu jest opracowanie i~implementacja metod związanych z detekcją i filtracją linii bazowej w sygnale EKG.
\item[Dane] Do przetwarzania w projekcie wykorzystywane są ciągi próbek sygnału EKG z bazy Physionet.org (\verb|MIT-BIH|)
\item[Szukane] Efektem pracy ma być moduł programu filtrujący sygnał EKG z zakłóceń sieciowych/mięśniowych oraz usuwający z sygnału falowanie linii izoelektrycznej przy wykorzystaniu różnych metod. Moduł powinien umożliwiać wybór metody detekcji i usuwania oraz ustawienie parametrów.
\end{description}

\subsubsection{Badania literaturowe}
\label{sec:baseline:papers}
Problem detekcji i usuwania linii bazowej jest szeroko omawiany zarówno w książkach związanych z tematyką przetwarzania sygnałów diagnostycznych \cite{AUGUST1}, jak i w publikacjach i pracach naukowych \cite{KAUR1}\cite{MUKA1}. W kontekście tego klasycznego problemu porównywano już wiele metod -- zarówno typowe filtry FIR i IIR, transformację falkową, uśrednianie sygnału, czy też filtry adaptacyjne, biorąc pod uwagę zmiany współczynników jakości w zależności od ustalonych parametrów. Wśród parametrów, jakie stosowane są do oceny jakości działania filtrów można wymienić:
\begin{enumerate}
\item \verb+SNR+  -- \textit{Signal-to-Noise Ratio} --  stosunek sygnału do szumu, jeden z popularniejszych parametrów
\item \verb+PSD+  -- \textit{Power Spectral Density} -- gęstość widmowa mocy
\item \verb+PRD+  -- \textit{Percent Root Mean Square Difference} -- procentowa różnica średnio-kwadratowa
\item \verb+PSNR+  -- \textit{Peak Signal-to-Noise Ratio} -- szczytowy stosunek sygnału do szumu
\end{enumerate}
Każda metoda ma zarówno wady i zalety, dlatego dokonanie wyboru w celu implementacji w projekcie nie było łatwym zadaniem.

\paragraph{Wybór metod}
\mbox{} \\
Biorąc pod uwagę fakt, iż nadrzędnym celem projektu jest zbudowanie oprogramowania do analizy sygnału EKG, a nie badanie wydajności metod na poszczególnych etapach przetwarzania, wybór padł na metody sprawdzone. Wybrano jedną metodę uśredniającą -- \verb|metodę średniej| \verb|kroczącej| (\textit{ang. moving average method}) oraz filtry o \verb|nieskończonej odpowiedzi impulsowej| (\verb|IIR|) -- \verb|Butterwortha| i \verb|Czebyszewa|. Taki zestaw, bazując na informacjach ze źródeł, powinien spełnić stawianie przed nim wymagania. Do oceny jakości działania filtrów zastosowano \verb|SNR| i \verb|PSNR|.

\subparagraph{Metoda średniej kroczącej}
\mbox{} \\
Algorytm uśredniania sygnału w ruchomym oknie jest bardzo prostą metodą usuwania zakłóceń sygnału. Sygnał jest wygładzany poprzez zastąpienie wartości średnią wyliczoną wartością sygnału w~oknie o ustalonej szerokości (nieparzystej, \textit{2N+1}, gdzie \verb|N| to sąsiedztwo) \eqref{eq:I}. Po dokonaniu obliczenia okno przesuwane jest o jedną próbkę i kalkulacje są powtarzane. 
\begin{equation}
Ys(i) = \frac{1}{2N+1}(Y(i+N)+Y(i+N-1)+...+Y(i-N) \label{eq:I}
\end{equation}
Oczywistym jest, że na krańcach przedziału nie można stosować tej metody, ze względu na brak możliwości obrania odpowiednio dużego sąsiedztwa -- w takim przypadku przepisuje się wartości sygnału bazowego do wektora wynikowego.\\
Aby średnią kroczącą zastosować do usuwania lini bazowej należy wyliczone wartości odejmować od odpowiednich próbek sygnału bazowego. Metoda zachowuje się wtedy jak swego rodzaju filtr górnoprzepustowy \cite{KAUR1}. Metoda ta została wybrana do projektu ze względu na prostotę i szybkość implementacji. Nie jest ona jednak pozbawiona wad. Usuwa zakłócenia w sposób zadowalający, niemniej jednak powoduje dość znaczące zatracanie cech charakterystycznych sygnału, na przykład w przypadku EKG wyraźnych "pików". Jest to łatwe do stwierdzenia empirycznie, zestawiając wykresy sygnału bazowego i przefiltrowanego \eqref{fig:fig_srednia}.

\subparagraph{Filtr Butterwortha}
\mbox{} \\
Według różnych źródeł \cite{JAGT1} pasmo informacyjne sygnału EKG zawiera się w przedziale od \verb|0.05Hz-5Hz| do \verb|100Hz|, dlatego do redukcji zakłóceń można stosować filtry zarówno dolno- jak i górnoprzepustowe.\\

Filtr \verb|Butterwortha| to filtr charakteryzujący się maksymalnie płaską charakterystyką amplitudową w paśmie przenoszenia. Częstotliwość graniczną filtru wyznacza spadek sygnału (\textit{ang. ripple}) o 3 dB. Jest to filtr powszechnie uważany za bardzo skuteczny w filtrowaniu linii bazowej \cite{KAUR1}.\\
Na potrzeby projektu proponowane są dwie wersjie filtru \verb|Butterwortha| (IIR):
\begin{enumerate}
\item \textit{górnoprzepustowa} drugiego rzędu, o częstotliwości odcięcia (\textit{ang. cut-off frequency}) \verb|5Hz|\eqref{fig_butter}
\item \textit{dolnoprzepustowa} czwartego rzędu, o częstotliwości odcięcia \verb|100Hz|
\end{enumerate}
O ile w przypadku filtru dolnoprzepustowego dobór częstotliwości nie dziwi -- bo jak już wcześniej wspomniano na 100Hz kończy się pasmo informacyjne EKG -- o tyle w przypadku górnoprzepustowego może paść pytanie, dlaczego akurat 5Hz. W sygnale EKG, do poziomu 5Hz, można napotkać artefakty wynikające z ruchy pacjenta (\textit{motion artifacts}, a dodatkowo filtr wyeliminuje również zakłócenia na niższych częstotliwościach.
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{include/fig_butter.eps}
\caption{Odpowiedź amplitudowa filtru Butterwortha (górnoprzepustowego)}
\label{fig_butter}
\end{figure}
Jak już to zostało wspomniane -- filtr Butterwortha jest bardzo skuteczny, niemniej jednak równniez on nie jest pozbawiony wad. Przez to, że jako filtr IIR ma nieliniową odpowiedź fazową, istnieje tendencja do znacznego zniekształcania sygnału bazowego. Aby uniknąć tego zjawiska stosuje się tzw. \verb|filtrowanie o charakterystyce zero-fazowej| \textit{ang. Zero-Phase filtering}, opisan w kolejnych paragrafch.

\subparagraph{Filtr Czebyszewa} 
\mbox{} \\
Jest to rodzaj filtru elektrycznego, którego charakterystyczną cechą jest wykorzystanie \textit{wielomianów Czebyszewa} do aproksymacji charakterystyki częstotliwościowej amplitudowej. Optymalizacja przebiegu charakterystyki częstotliwościowej amplitudowej w filtrach tego typu ma kluczowe znaczenie.\\
Wyróżnia się dwa typy filtrów:
\begin{enumerate}
\item \verb|filtry Czebyszewa typu I| -- zafalowania przebiegu wzmocnienia w paśmie przepustowym, oraz płaski przebieg charakterystyki w paśmie zaporowym
\item \verb|filtry Czebyszewa typu II| -- zafalowania przebiegu wzmocnienia w paśmie zaporowym, oraz płaski przebieg charakterystyki w paśmie przepustowym, nazywany przez to zwykle \textit{inwersyjnym}
\end{enumerate}
Filtry \verb|Czebyszewa| mają porównywalną skuteczność do filtrów \verb|Butterwortha|. Na ich korzyśc przemawia zdecydowanie lepsze tłumienie za pasmem przepustowym (stroma charakterystyka). Wsród minusów można wymienić falowanie w pasmie przepustowym w przypadku filtru \verb|typu I| oraz dodatkowy parametru -- \verb|tętnienie, ang. ripple| -- niezbędny do zdefiniowania filtru. Dla porównania w~przypadku \verb|Butterwortha| wystarczy częstotliwośc odcięcia.\\
Na potrzeby projektu proponowane są dwie wersjie filtru \verb|Czebyszewa|:
\begin{enumerate}
\item \textit{górnoprzepustowa typu II}, trzeciego rzędu, o częstotliwości odcięcia \verb|0.5Hz| i tętnieniu \verb|-3dB|
\item \textit{dolnoprzepustowa typu I}, piątego rzędu, o częstotliwości odcięcia \verb|100Hz| i tętnieniu \verb|-3dB|
\end{enumerate}
Również w przypadku tego filtru konieczne będzie zastosowanie się tzw. \verb|filtrowania o charakterystyce| \verb|zero-fazowej| \textit{ang. Zero-Phase filtering}.

\subparagraph{Filtrowanie o charakterystyce zero-fazowej}
\mbox{} \\
Filtrowanie \verb|Zero-Phase| wykorzystuje w swym działaniu współczynniki filtru bazowego (np. Butterwortha czy też Czebyszewa) i wykonuje filtrowanie dwukrotnie, stosując małe okno, w obie strony sygnału -- od początku do końca i od końca do początku. Rezultatem tego filtrowania jest:
\begin{enumerate}
\item brak opóźnień w sygnale przefiltrowanym, ze względu na przejście sygnału w obie strony
\item redukcja zaburzeń sygnału na krańcach przedziału -- poprzez zachowanie warunków początkowych
\item zachowanie cech charakterystycznych sygnału bazowego (w przeciwieństwie do samego filtru IIR czy też np. średniej kroczącej)
\item rząd filtru zero-fazowego jest równy podwojonemu rzędowi filtru bazowego.
\end{enumerate}
Filtr może być opisany przez poniższe rekurencyjne równanie \eqref{eq:II}
\begin{equation}
y[n] = a_{0}x[n]+a_{1}x[n+1]+a_{2}x[n+2]+a_{3}x[n+3]+b_{1}y[n+1]+b_{2}y[n+2]+b_{3}y[n+3] \label{eq:II}
\end{equation}
a nastepnie odwrotne rekurencyjne równanie, do zastosowania przy przejściu od końca sygnału \eqref{eq:III}
\begin{equation}
y[n] = a_{0}x[n]+a_{1}x[n-1]+a_{2}x[n-2]+a_{3}x[n-3]+b_{1}y[n-1]+b_{2}y[n-2]+b_{3}y[n-3] \label{eq:III}
\end{equation}
Należy przyjąć, że \verb|b| to współczynniki licznika transmitancji filtru bazowego, a \verb|a| to współczynniki mianownika.

\subparagraph{Signal-to-Noise Ratio}
\mbox{} \\
Parametr oznaczający stosunek sygnału (użytecznego) do szumu. Jest stosowany jako parametr pomocniczy w określeniu jakości sygnału. Obliczany za pomocą następującej formuły:
\begin{equation}
SNR= 20\log_{10}(\frac{A_{sygnal}}{A_{szum}}) \label{eq:IV}
\end{equation}
gdzie \verb|A| oznacza amplitudę sygnału. W praktyce \verb|A| wyliczane jest jako średnia kwadratowa amplitud w sygnale.

\subparagraph{Peak Signal-to-Noise Ratio}
\mbox{} \\
Parametr oznaczający stosunek maksymalnej mocy sygnału do średniej wartości szumu. Również jest stosowany jako parametr pomocniczy w określeniu jakości sygnału. Obliczany za pomocą następującej formuły:
\begin{equation}
PSNR= 10\log_{10}(\frac{k^2}{MSE}) \label{eq:V}
\end{equation}
gdzie \verb|MSE| oznacza błąd średniokwadratowy obliczany za pomoca
\begin{equation}
MSE=\frac{1}{N}\sum\limits_{i=1}^N([f(i)-f'(i)]^2)\label{eq:VI}
\end{equation}
gdzie \verb|N| oznacza długość sygnału, \verb|f| to sygnał oryginalny, \verb|f'| to sygnał przefiltrowany, a \verb|k| to maksymalna amplituda sygnału bazowego.


\subsubsection{Opis funkcjonalności modułu}
Nadrzędnym przeznaczeniem modułu jest filtracja i usuwanie lini bazowej. Niemniej jednak moduł bardziej kompleksowo przygotowuje sygnał do dalszej analizy.\\
Realizowane są następujące funkcjonalności:
\begin{enumerate}
\item Wybór kanału, który będzie przetwarzany w kolejnych etapach działania aplikacji. Wybór dokonywany jest poprzez porównanie SNR obu kanałów.
\item Detekcja i usunięcie zakłóceń lini bazowej sygnału -- główne zadanie.
\item Normalizacja sygnału przed zwróceniem do dalszej analizy.
\end{enumerate}
\begin{figure}[h]
\centering
\includegraphics[scale=0.8]{include/baselineflow.eps}
\caption{Diagram prezentujący funkcjonalnosć modułu ECGBaselineRemoval}
\label{fig:fig_baselineflow}
\end{figure}
Koniecznosć wprowadzenia oceny i wyboru kanału do dalszej analizy został podyktowany nadrzędnymi decyzjami architekturalnymi. W razie potrzeby może zostać pominięta.\\

\subsubsection{Opis procedur i metod}
\label{sec:baseline:procs}
\paragraph{Klasa: ECGBaselineRemoval}
- Główna klasa modułu, dziedzicząca po wirtualnym ECGBaselineModule. Jest ona odpowiedzialna za przetworzenie sygnału pierwotnego i zwrócenie sygnału przefiltrowanego (pozbawionego zakłóceń linii bazowej). Pojedyncza instancja tej klasy tworzona jest wraz z uruchomieniem aplikacji.

\begin{lstlisting}
void ECGBaselineRemoval::runModule(const ECGSignal &inputSignal, 
const ECGInfo & ecgi, ECGSignalChannel &outputSignal)
\end{lstlisting}
Podstawowa metoda w module. Jest odpowiedzialna za uruchomienie głównej funkcjonalności -- przetworzenia sygnału. Pod hasłem "przetworzenie" kryje się wybranie odpowiedniego kanału, realizacja filtrowania oraz normalizacja. \\
Parametry:
\begin{enumerate}
\item \verb|const ECGSignal &inputSignal| -- referencja do sygnału bazowego, zawierającego dwa kanały
\item \verb|const ECGInfo & ecgi| -- referencja do obiektu informacyjnego typu \verb|ECGInfo|, przechowującego dane dotyczące sygnału (m.in. częstotliwość próbkowania, ilość próbek)
\item \verb|ECGSignalChannel &outputSignal| -- referencja do sygnału wyjściowego -- po wykonaniu metody \verb|ECGBaselineRemoval::runModule| odnosi się do sygnału przefiltrowanego; instancja klasy \verb|ECGSignalChannel| -- pojedynczy kanał
\end{enumerate}

\begin{lstlisting}
void ECGBaselineRemoval::setParams(ParametersTypes &params)
\end{lstlisting}
Metoda pozwala na wybór algorytmu filtracji poprzez ustawienie flagi \verb|baselineRemovalMethod|, oraz ustawienie parametrów, korzystając z wartości pochodzących z graficznego interfejsu użytkownika. \\
Parametry:
\begin{enumerate}
\item \verb|ParametersTypes &params| -- referencja do obiektu typu \verb|PatametersTypes| -- zawierającego parametry przekazane z graficznego interfejsu użytkownika
\end{enumerate}

\begin{lstlisting}
void ECGBaselineRemoval::movingAverageBaselineRemoval(ECGSignalChannel 
&inputSignal, ECGSignalChannel &outputSignal, int span)
\end{lstlisting}
Metoda opakowująca (tzw. \textit{wrapper}) instancję klasy realizującej filtrowanie metodą średniej kroczącej \ref{sec:baseline:papers}. \\
Parametry parametry:
\begin{enumerate}
\item \verb|ECGSignalChannel &inputSignal| -- referencja do obiektu będącego wybranym kanałem z sygnału wejściowego
\item \verb|ECGSignalChannel &outputSignal| -- referencja do obiektu wyjściowego
\item \verb|int span| -- wartość całkowita -- nieparzysta -- będąca rozmiarem okna w algorytmie średniej kroczącej
\end{enumerate}

\begin{lstlisting}
void ECGBaselineRemoval::butterworthBaselineRemoval(ECGSignalChannel 
&inputSignal, ECGSignalChannel &outputSignal, const ECGInfo &ecgInfo, 
int order, double cutoffFrequency, double attenuation)
\end{lstlisting}
Metoda opakowująca instancję klasy realizującej filtrowanie filtrem Butterwortha \ref{sec:baseline:papers}. \\

Parametry:
\begin{enumerate}
\item \verb|ECGSignalChannel &inputSignal| -- referencja do obiektu będącego wybranym kanałem z sygnału wejściowego
\item \verb|ECGSignalChannel &outputSignal| -- referencja do obiektu wyjściowego
\item \verb|const ECGInfo &ecgi| -- referencja do obiektu informacyjnego \verb|ECGInfo|, przechowującego dane sygnału (m.in. częstotliwość próbkowania, ilość próbek)
\item \verb|int order| -- wartość całkowita dodatnia, odpowiadająca rzędowi filtra
\item \verb|double cutoffFrequency| -- częstotliwość odcięcia filtra -- wartość zmiennoprzecinkowa (może przyjąć np \textit{0.05[Hz]})
\item \verb|double attenuation| -- tłumienie 
\end{enumerate}

\begin{lstlisting}
void ECGBaselineRemoval::chebyshevBaselineRemoval(ECGSignalChannel 
&inputSignal, ECGSignalChannel &outputSignal, const ECGInfo &ecgInfo, 
int order, double cutoffFrequency, double ripple)
\end{lstlisting}
Metoda opakowująca instancję klasy realizującej filtrowanie filtrem Czebyszewa \ref{sec:baseline:papers}. \\

Parametry:
\begin{enumerate}
\item \verb|ECGSignalChannel &inputSignal| -- referencja do obiektu będącego wybranym kanałem z sygnału wejściowego
\item \verb|ECGSignalChannel &outputSignal| -- referencja do obiektu wyjściowego
\item \verb|const ECGInfo &ecgi|  -- referencja do obiektu informacyjnego \verb|ECGInfo|, przechowującego dane sygnału (m.in. częstotliwość próbkowania, ilość próbek)
\item \verb|int order|  -- wartość całkowita dodatnia, odpowiadająca rzędowi filtra
\item \verb|double cutoffFrequency|  -- częstotliwość odcięcia filtra  -- wartość zmiennoprzecinkowa (może przyjąć np \textit{0.05[Hz]})
\item \verb|double ripple|  -- tętnienie 
\end{enumerate}

\paragraph{Klasa: SignalPreprocessor}
- Klasa pomocnicza odpowiedzialna za wstępną analizę sygnału oraz jego normalizację.

\begin{lstlisting}
int SignalPreprocessor::evaluateSignalChannels(const ECGSignal &inputSignal)
\end{lstlisting}
Metoda realizująca porównanie kanałów, bazując na wyliczeniu współczynnika \verb|SNR| (\verb|Snr::snr|) tysiąca pierwszych próbek oraz wybór lepszego z nich. \\

Parametry:
\begin{enumerate}
\item \verb|const ECGSignal &inputSignal|  -- referencja do sygnału bazowego (dwa kanały)
\end{enumerate}
Wartość zwracana:
\begin{enumerate}
\item liczba całkowita odpowiadająca numerowi kanału po ewaluacji (1 albo 2)
\end{enumerate} 

\begin{lstlisting}
void SignalPreprocessor::normalize(ECGSignalChannel &inputSignal, 
ECGSignalChannel &outputSignal)
\end{lstlisting}
Metoda odpowiedzialna za normalizację sygnału, tj. przeskalowanie do przedziału wartości 0..1. Wywoływana jako ostatni etap realizacji przetwarzania wstępnego  -- \verb|ECGBaselineRemoval::runModule()|.\\
Parametry:
\begin{enumerate}
\item \verb|ECGSignalChannel &inputSignal|  -- referencja do obiektu wejściowego
\item \verb|ECGSignalChannel &outputSignal|  -- referencja do obiektu wyjściowego
\end{enumerate}

\paragraph{Klasa: Filter}
- Klasa odpowiadająca za realizację filtrowania \verb|IIR| (\textit{ang. infinite impulse response}) na bazie współczynników pozyskanych z metod projektujących filtry \verb|Butterwortha| bądź \verb|Czebyszewa|. 

\begin{lstlisting}
void Filter::zeroPhase(std::vector<double> b, std::vector<double> a,
ECGSignalChannel &inputSignal, ECGSignalChannel &outputSignal, int order)
\end{lstlisting}
Metoda realizująca filtrowanie typu Zero-Phase \ref{sec:baseline:papers}. Wykorzystuje współczynniki transmitancji filtru bazowego, zgodnie z \eqref{eq:II} oraz \eqref{eq:III}. Metoda została zaimplementowana na wzór funkcji \verb|filtfilt| ze środowiska Matlab.\\

Parametry:
\begin{enumerate}
\item \verb|(std::vector<double> b|  -- wektor współczynników pochodzących z licznika transmitancji zaprojektowanego filtra
\item \verb|(std::vector<double> a|  -- wektor współczynników pochodzących z mianownika transmitancji zaprojektowanego filtra
\item \verb|ECGSignalChannel &inputSignal|  -- referencja do obiektu wejściowego
\item \verb|ECGSignalChannel &outputSignal|  -- referencja do obiektu wyjściowego
\item \verb|int order|  -- wartość całkowita dodatnia, odpowiadająca rzędowi filtra
\end{enumerate}

\begin{lstlisting}
void Filter::filter(std::vector<double> b, std::vector<double> a,
ECGSignalChannel &inputSignal, ECGSignalChannel &outputSignal, int order)
\end{lstlisting}

Metoda realizująca filtr \verb|IIR| na bazie zestawu współczynników transmitancji filtru \verb|Butterwortha| bądź \verb|Czebyszewa|. Funkcja bazuje na swoim odpowiedniku (\verb|filter|) ze środowiska Matlab. Jest uniwersalna i może realizować również inne filtry. \verb|Filter::filter| może być wykorzystywana autonomicznie bądź wywoływana niejawnie w ramach metody \verb|Filter::zeroPhase|.\\

Parametry:
\begin{enumerate}
\item \verb|(std::vector<double> b|  -- wektor współczynników pochodzących z licznika transmitancji zaprojektowanego filtra
\item \verb|(std::vector<double> a|  -- wektor współczynników pochodzących z mianownika transmitancji zaprojektowanego filtra
\item \verb|ECGSignalChannel &inputSignal|  -- referencja do obiektu wejściowego
\item \verb|ECGSignalChannel &outputSignal|  -- referencja do obiektu wyjściowego
\item \verb|int order|  -- wartość całkowita dodatnia, odpowiadająca rzędowi filtra
\end{enumerate}

\paragraph{Klasa: MovingAverage}
- Klasa realizująca filtrowanie metodą średniej kroczącej.

\begin{lstlisting}
void MovingAverage::removeBaseline(ECGSignalChannel &inputSignal, 
ECGSignalChannel \&outputSignal, int span)
\end{lstlisting}
Metoda realizująca uswanie linii bazowej algorytmem średniej kroczącej. Główna pętla metody implementuje \eqref{eq:I} wywołując metodę \verb|MovingAverage::calculateAvgValueOfNeighbours|\\

Parametry:
\begin{enumerate}
\item \verb|ECGSignalChannel &inputSignal|  -- referencja do obiektu będącego wybranym kanałem z sygnału wejściowego
\item \verb|ECGSignalChannel &outputSignal|  -- referencja do obiektu wyjściowego
\item \verb|int span|  -- wartość całkowita  -- nieparzysta  -- będąca rozmiarem okna w algorytmie średniej kroczącej
\end{enumerate}

\begin{lstlisting}
double MovingAverage::calculateAvgValueOfNeighbours(gsl_vector *signal, 
int currentIndex, int span)
\end{lstlisting}
Metoda obliczająca średnią wartość próbek w sąsiedztwie danej próbki, w oknie o szerokości określonej wartością zmiennej \verb|int span|\\

Przyjmowane parametry
\begin{enumerate}
\item \verb|gsl_vector *signal|  -- wektor przechowujący próbki przetwarzanego sygnału
\item \verb|int currentIndex|  -- indeks aktualnie analizowanej próbki
\item \verb|int span|  -- wartość całkowita  -- nieparzysta  -- będąca rozmiarem okna w algorytmie średniej kroczącej
\end{enumerate}

\paragraph{Klasa: Butterworth}
- Klasa realizująca projektowanie filtra \verb|Butterwortha|.

\begin{lstlisting}
std::vector<std::vector<double>> Butterworth::filterDesign(int order, 
double attenuation, double cutoffFreq, int sampleFreq, int type)
\end{lstlisting}
Metoda odpowiedzialna za zaprojektowanie filtra \verb|Butterwortha|.\\

Parametry:
\begin{enumerate}
\item \verb|int order|  -- wartość całkowita dodatnia, odpowiadająca rzędowi filtra
\item \verb|double attenuation|  -- tłumienie 
\item \verb|double cutoffFrequency|  -- częstotliwość odcięcia filtra  -- wartość zmiennoprzecinkowa (może przyjąć np \textit{0.05[Hz]})
\item \verb|int sampleFreq|  -- częstotliwość próbkowania sygnału wejściowego, potrzebną do obliczenia współczynników filtra
\item \verb|int type|  -- zmienna określająca typ (lowpass -1, highpass - 2) filtra do zaprojektowania
\end{enumerate}
Wartość zwracana:
\begin{enumerate}
\item \verb|std::vector<std::vector<double>>|  -- dwuwymiarowy wektor wartości zmiennoprzecinkowych, zwracający zestaw współczynników z mianownika i licznika transmitancji zaprojektowanego filtra
\end{enumerate}

\paragraph{Klasa: Chebyshev}
- Klasa realizująca projektowanie filtra \verb|Czebyszewa|.\\

\begin{lstlisting}
std::vector<std::vector<double>> Chebyshev::filterDesign(int order, 
double ripple, double cutoffFreq, int sampleFreq, int type)
\end{lstlisting}
Metoda odpowiedzialna za zaprojektowanie filtra \verb|Czebyszewa|\\

Parametry:
\begin{enumerate}
\item \verb|int order|  -- wartość całkowita dodatnia, odpowiadająca rzędowi filtra
\item \verb|double ripple|  -- tętnienie pasma przepustowego 
\item \verb|double cutoffFrequency|  -- częstotliwość odcięcia filtra  -- wartość zmiennoprzecinkowa (może przyjąć np \textit{0.05[Hz]})
\item \verb|int sampleFreq|  -- częstotliwość próbkowania sygnału wejściowego, potrzebną do obliczenia współczynników filtra
\item \verb|int type|  -- zmienna określająca typ (\textit{lowpass} - 1, \textit{highpass} - 2) filtra do zaprojektowania
\end{enumerate}
Wartość zwracana:
\begin{enumerate}
\item \verb|std::vector<std::vector<double>>|  -- dwuwymiarowy wektor wartości zmiennoprzecinkowych, zwracający zestaw współczynników z mianownika i licznika transmitancji zaprojektowanego filtra
\end{enumerate}

\paragraph{Klasa: Snr}
- Klasa realizująca metody ewaluacji filtrów, wraz z metodami pomocniczymi zgodnie ze wzorami \eqref{eq:IV},\eqref{eq:V} i \eqref{eq:VI}. Pozwala na wyliczenie współczynników \verb|SNR| (\verb|Snr::snr|) oraz \verb|PSNR| (\verb|Snr::psnr|) dla zadanych sygnałów. Ze względu na generyczność rozwiązań w niej zaimplementowanych opis interfejsu zostaje pominięty.

\paragraph{Projektowanie filtrów  -- realizacja}\mbox{} \\
Do projektowania filtrów \verb|Butterwortha| oraz \verb|Czebyszewa| miała zostać wykorzystana biblioteka \verb|SPUC -| \verb|Signal Processing using C++| (\textit{http://spuc.sourceforge.net/}). Niestety po jej zbudowaniu zgodnie z zaleceniami twórców okazało się, że praktycznie niemożliwe jest linkowanie jej w sposób klasyczny z~projektem Visual Studio. Ze względu na ten fakt próby jej wykorzystania zostały porzucone na rzecz implementacji analogicznej z tą, która występuje w pakiecie Matlab.

\subsubsection{Warunki testowania}
\label{sec:baseline:tests}
Testy na etapie prototypowania przeprowadzane były w środowisku Matlab. Dało to pogląd na wymagania, jakie niesie za sobą proces usuwania linii bazowej. Na etapie realizacji projektu przeprowadzano równolegle testy zarówno w aplikacji jak i ponownie w Matlab. Dane testowe pochodzą z~bazy \verb|MIT-BIH|.\\\\
Wartości proponowanych parametrów dla poszczególnych metod zostały dobrane na bazie doświadczeń oraz informacji ze źródeł naukowych, o czym zostało już nadmienione powyżej \ref{sec:baseline:papers}.\\

\paragraph{Parametry filtrów}
\subparagraph{Metoda średniej kroczącej}
\begin{enumerate}
\item \verb|span = 5|
\end{enumerate}

\subparagraph{Filtr Butteworth'a}
\mbox{} \\
Pierwszy zestaw:
\begin{enumerate}
\item \verb|order = 2|
\item \verb|cutoff_frequency = 5Hz|
\item \verb|sample_frequency = 360Hz|
\item \verb|type = 'high'|
\end{enumerate}
Drugi zestaw:
\begin{enumerate}
\item \verb|order = 4|
\item \verb|cutoff_frequency = 100Hz|
\item \verb|sample_frequency = 360Hz|
\item \verb|type = 'low'|
\end{enumerate}

\subparagraph{Filtr Czebyszewa}
\mbox{} \\
Pierwszy zestaw:
\begin{enumerate}
\item \verb|order = 3|
\item \verb|cutoff_frequency = 0.5Hz|
\item \verb|ripple = 3dB|
\item \verb|sample_frequency = 360Hz|
\item \verb|type = 'high'|
\end{enumerate}
Drugi zestaw:
\begin{enumerate}
\item \verb|order = 5|
\item \verb|cutoff_frequency = 100Hz|
\item \verb|ripple = 3dB|
\item \verb|sample_frequency = 360Hz|
\item \verb|type = 'low'|
\end{enumerate}

W większości przypadków wyodrębnione w procesie prototypowania skuteczne parametry pokryły się z propozycjami naukowców.

\subsubsection{Wyniki}
\label{sec:baseline:results}
Testy polegały na analizie poglądowej zmiany przebiegu sygnału po zastosowaniu filtra oraz na wyliczeniu parametrów \verb|SNR| oraz \verb|PSNR|.

\paragraph{Poglądowe porównanie metod}
\subparagraph{Metoda średniej kroczącej}
\mbox{} \\
W przypadku jedynej metody uśredniającej w zestawieniu wynik poglądowy jest zadowalający. Na wykresie \ref{fig:fig_srednia} widać, iż samo zastosowanie uśredniania wygładziło przebieg sygnału (wykres pomarańczowy), a po odjęciu wartości wyliczonej od sygnału bazowego otrzymano sygnał wynikowy wydaje się być bardziej podatny na skuteczną analizę.\\

\begin{figure}[h]
\centering
\includegraphics[scale=0.8]{include/fig_srednia.eps}
\caption{Sygnał bazowy i sygnał przefiltrowany algorytmem średniej kroczącej}
\label{fig:fig_srednia}
\end{figure}

Niestety wyraźnie widać utratę cech charaktetystycznych sygnału bazowego  -- "piki" dalej występują w tych samych miejscach, niemniej jednak zostały one wyraźnie spłaszczone. Istnieje podejrzenie, że nie każda metoda ekstrakcji cech sygnału poradzi sobie z tak niewyraźnymi "pikami". W tym przypadku należy przeprowadzić testy polegające na stosowaniu różnych metod analizy sygnału, aby ostatecznie ocenić skuteczność średniej kroczącej.

\subparagraph{Filtry górnoprzepustowe}
\mbox{} \\
Filtry górnoprzepustowe, niezależnie od typu -- czy \verb|Butterwortha|, czy też \verb|Czebyszewa|  -- wykazały się dobrą skutecznoscią w usuwaniu linii bazowej. Jak widać na wykresie \ref{fig:fig_skuteczne} falowanie sygnału zostało w~znaczącym stopniu zredukowane. Minimalnie lepszy efekt w przypadku filtrów górnoprzepustowych dawał \verb|Butterworth|.

\begin{figure}[h]
\centering
\includegraphics[scale=0.8]{include/fig_skuteczne.eps}
\caption{Sygnał bazowy i sygnały po skutecznej filtracji}
\label{fig:fig_skuteczne}
\end{figure}

Warto nadmienić, że filtry \verb|IIR| nie powodują aż tak dużej utraty cech charakterystycznych jak średnia krocząca. 

\subparagraph{Filtry dolnoprzepustowe}
\mbox{} \\
Bez względu na to, że wartości filtrów dolnoprzepustowych były bezpośrednio zapożyczone z literatury, ich stosowanie nie przyniosło oczekiwanych rezultatów. Przebieg sygnału po przefiltrowaniu pozostawał praktycznie bez zmian. Jedyną korzyścią, jaką można dostrzec jest nieznaczne wygładzenie wykresu (\ref{fig:fig_nieskuteczne}).

\begin{figure}[h]
\centering
\includegraphics[scale=0.8]{include/fig_nieskuteczne.eps}
\caption{Sygnał bazowy i sygnały po nieskutecznej filtracji}
\label{fig:fig_nieskuteczne}
\end{figure}

\paragraph{Porównanie parametrów SNR i PSNR}
\mbox{} \\
Oprócz oceny poglądowej działania filtrów dokonano również pomiarów parametrów \verb|SNR| i \verb|PSNR|. Zestawienie wyników prezentują poniższe tabele:

\begin{center}
\verb|SNR|[dB]\\ 

\begin{tabular}{|c|c|c|c|c|c|}
\hline
sygnał & Butter. (high) & Czeby. (high) & Śr. krocząca & Butter. (low) & Czeby. (low) \\ \hline
100\_V5   & -6.4236 & -3.9277 & -12.5362  & 36.7097 &  18.4894 \\ \hline
100\_MII   & -6.0413 &  -4.2806 &  -12.7211 & 40.6467 &  18.1320 \\ \hline
\end{tabular}
\mbox{} \\
\mbox{} \\
\verb|PSNR|[dB]\\
\begin{tabular}{|c|c|c|c|c|c|}
\hline
sygnał & Butter. (high) & Czeby. (high) & Śr. krocząca & Butter. (low) & Czeby. (low) \\ \hline
100\_V5   & -48.1567 & -61.0092 & -41.1130  &  -75.2726 &  -56.3931 \\ \hline
100\_MII   &  -50.9663 &   -76.0654 &  -43.5110 &  -72.6385 &  -57.2145 \\ \hline
\end{tabular}
\mbox{} \\
\end{center}
Niestety oba parametry w przypadku badanych sygnałów okazały się zwracać nieoczekiwane  -- pozornie błędne  -- wyniki. W idealnym przypadku zarówno \verb|SNR| jak i \verb|PSNR| powinny zwracać wartości dodatnie. Co więcej  -- bardziej skuteczne metody (w tym przypadku filtry górnoprzepustowe oraz metoda średniej kroczącej) powinny zwracać wyniki znacząco wyższe od metod nieskutecznych. Otrzymane wyniki mogłyby wskazywać na istnienie błędu w implementacji metody, jednak była ona weryfikowana wielokrotnie.\\
Dopóki nie zostaną roziwązane zaistniałe problemy nie można traktować parametrów \verb|SNR| oraz \verb|PSNR| jako wiążących w procesie oceny skuteczności filtracji.

\paragraph{Wnioski}
\mbox{} \\
Zaproponowane metody w przypadku doboru odpowiednich parametrów oferują skuteczność na dobrym poziomie. Potwierdzają to również lepsze rezultaty osiągane przez zespoły odpowiedzialne za moduły analizujące sygnał EKG po zmianie sygnału źrodłowego na sygnał pozbawiony lini bazowej.\\

Testy wykazały, że filtry górnoprzepustowe mogą z powodzeniem być stosowane w tym module. Metoda średniej kroczącej do pewnego stopnia również dobrze poradziła sobie z zadaniem. Warto jednak zaznaczyć, że w przypadku wyboru między tymi dwoma typami filtracji i usuwania lini bazowej użytkownik musi zastanowić się, czy ważniejsze w dalszym etapie analizy będzie zachowanie wyraźnych cech sygnału czy też jak najskuteczniejsza redukcja falowania linii izoelektrycznej.\\

Zalecane są dalsze testy, w celu ustalenia przyczyn niskiej skuteczności filtrów dolnoprzepustowych.


\subsection{Wykrywanie załamków R}
\label{sec:Rs}
Autorzy: Paweł Maślanka i Norbert Pabian

\subsubsection{Opis zadania}
\label{sec:rs:desc}
\begin{description}
\item[Temat] Detekcja załamka R
\item[Opis] Detekcja załamków R w przefiltrowanym sygnale cyfrowym EKG, mimo iż jest tematem przebadanym i dobrze znanym, to wciąż nie istnieje złoty standard postępowania. Alternatywą dla powszechnie znanego algorytmu "Pan-Tompkins", opartego o podstawowe przekształcenia matematyczne potęgowania i różniczkowania może okazać się transformata Hilberta odpowiednio przygotowanego sygnału. Celem projektu jest opracowanie i implementacja metod służących do detekcji załamków R w sygnale EKG w oparciu o algorytm "Pan-Tompkins" oraz zespolone przekształcenie Hilberta. Podejmując projekt należy zwrócić uwagę na problematykę wystąpienia szumu i analizy sygnału na końcach przedziału, gdzie rozważania teoretyczne transformaty Hilberta nieco odbiegają od praktycznych.
\item[Dane] przefiltrowany ciąg próbek z modułu ECG\_BASELINE
\item[Szukane] moduł programu wyznaczający numery próbek załamków R oraz zaznaczający na
wykresie wykryte załamki R. Numery próbek pozwolą na dalszą analizę ilościową
i jakościową sygnału EKG.
\end{description}

\subsubsection{Badania literaturowe}
\label{sec:rs:papers}


\paragraph{Usuwanie stałej DC i normalizacja sygnału}
Często w sygnałach EKG można spotkać dodatkową stałą DC w sygnale którą można skutecznie usunąć stosując wzór
\eqref{eq:a}. Obliczanie średniej sygnału i odejmowanie jej od każdej próbki sygnału jest przybliżoną metodą usuwania składowej stałej.
\begin{equation}
y(n) = x(n)-\frac{x(1)+x(2)+...+x(n)}{n} \label{eq:a}
\end{equation}
\indent Aby ułatwić wszystkie obliczenia dokonywane na wybranym sygnale można znormalizować sygnał tak by wartości amplitudy były w zakresie od $-1$ do $1$. Normalizację sygnału można uzyskać stosując wzór \eqref{eq:b}.
\begin{equation}
y(n) = \frac{x(n)}{\max(|x(n)|)} \label{eq:b} 
\end{equation}
\indent Rezultaty operacji usuwania składowej DC oraz normalizacji sygnału możemy zobaczyć na rysunku \eqref{znorm}.
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{include/surowy.eps}
\caption{Sygnał EKG oryginalny(100\_ MLII.dat)}
\label{surowy}
\end{figure}
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{include/znormalizowany.eps}
\caption{Sygnał EKG znormalizowany bez składowej stałej(100\_ MLII.dat)}
\label{znorm}
\end{figure}
\paragraph{Filtr dolnoprzepustowy}
Filtracja dolnoprzepustowa $0Hz - 15Hz$ ma na celu usunięcie zakłóceń sygnału powstających
pod wpływem drżenia mięśni ($35 Hz$), wpływu
zakłóceń sieci elektroenergetycznej $50/60 Hz$, wpływu załamka $T$ oraz pływającej izolinii
elektrycznej. Powszechnie stosowany jest filtr o transmitancji \eqref{eq:c}.
\begin{equation}
H(z) = \frac{1}{32} \frac{(1-z^{-6})^2}{(1-z^1)^2}\label{eq:c}
\end{equation}
\indent Dla częstotliwości próbkowania $fs = 200 Hz$, częstotliwość odcięcia
 dla filtru \eqref{eq:c} wynosi $fc = 11 Hz$. Filtr wprowadza opóźnienie $6$ próbek ($30 ms$). Dla $60 Hz$ filtr
tłumi sygnał na poziomie większym niż $35 dB$.
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{include/lowpassampres.eps}
\caption{Odpowiedź amplitudowa filtru dolnoprzepustowego}
\label{lowpassampres}
\end{figure}
\begin{figure}[htl]
\centering
\includegraphics[scale=0.5]{include/lowpassphaseres.eps}
\caption{Odpowiedź fazowa filtru dolnoprzepustowego}
\label{lowpassphaseres}
\end{figure}
\newpage
Rysunek \eqref{lowpassampres} oraz \eqref{lowpassphaseres} przedstawiają dokładne charakterystyki filtru dolnoprzepustowego. Filtr ten ma liniową odpowiedź fazową. Bardzo dobrze tłumi zakłócenia z linii energetycznych($50/60Hz$) oraz pozostałe zakłócenia na poziomie powyżej $25dB$.\\
\indent Rysunek \eqref{lowpass} przedstawia sygnał EKG po filtracji filtrem dolnoprzepustowym. Rezultat można porównać z oryginalnym sygnałem znajdującym się na rysunku \eqref{surowy}, możemy dostrzec znaczną eliminację zakłóceń pochodzących od wyższych częstotliwości.
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{include/lowpass.eps}
\caption{Sygnał EKG po filtracji dolnoprzepustowej (100\_ MLII.dat)}
\label{lowpass}
\end{figure}\\
\paragraph{Filtr górnoprzepustowy}
Filtracja górnoprzepustowa ma na celu wyeliminowanie wszystkich niepożądanych sygnałów o niskich częstotliwościach. Transmitancja użytego filtru jest przedstawiona równaniem \eqref{eq:d}
\begin{equation}
H(z) = 32z^{-16} -\frac{1-z^{-32}}{1-z^{-1}}\label{eq:d}
\end{equation}
\indent Dla częstotliwości próbkowania $fs = 200 Hz$, częstotliwość odcięcia wynosi $fc = 5 Hz$.
Filtr wprowadza opóźnienie $16$ próbek ($80 ms$).
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{include/hipassampres.eps}
\caption{Odpowiedź amplitudowa filtru górnoprzepustowego}
\label{hipassampres}
\end{figure}
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{include/hipassphaseres.eps}
\caption{Odpowiedź fazowa filtru górnoprzepustowego}
\label{hipassphaseres}
\end{figure}
\newpage
Rysunek \eqref{hipassampres} oraz \eqref{hipassphaseres} przedstawiają dokładne charakterystyki filtru górnoprzepustowego. 
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{include/hipass.eps}
\caption{Sygnał EKG po filtracji górnoprzepustowej(100\_ MLII.dat)}
\label{hipass}
\end{figure}\\
\paragraph{Różniczkowanie sygnału}
Po filtracji sygnału EKG, następnym krokiem jest różniczkowanie sygnału. Efektem takiego różniczkowania jest podkreślenie w sygnale nachyleń zespołów $QRS$ oraz stłumienie załamków $P$ i $T$. Przykładowa transmitancja do wykonania tej operacji jest przedstawiona równaniem \eqref{eq:e}.
\begin{equation}
H(z) = 0.1(2+z^{-1} - z^{-3} - 2z^{-4})\label{eq:e}
\end{equation}
\indent Na rysunku \eqref{dif} można dostrzec efekt różniczkowania sygnału.
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{include/diferentiation.eps}
\caption{Sygnał EKG po różniczkowaniu (100\_ MLII.dat)}
\label{dif}
\end{figure}\\
\paragraph{Potęgowanie sygnału}
Potęgowanie sygnału ma na celu jeszcze bardziej stłumić załamki $P$ i $T$ oraz wzmocnić próbki
sygnału reprezentujące zespół QRS. Dodatkowo odwraca ujemną część sygnału i dzięki temu ujemne zespoły $QRS$ również zostaną poddane analizie.
\begin{equation}
y(n) = (x(n))^2\label{eq:f}
\end{equation}
\indent Efekt potęgowania sygnału możemy zaobserwować na rysunku \eqref{squaring}.
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{include/squaring.eps}
\caption{Sygnał EKG po potęgowaniu (100\_ MLII.dat)}
\label{squaring}
\end{figure}\\
\newpage
\paragraph{Całkowanie sygnału w ruchomym oknie}
Całkowanie sygnału w ruchomym oknie ma na celu uzyskanie pojedynczej “fali” w obrębie
zespołu $QRS$. Bardzo ważną sprawą jest poprawne dobranie długości okna całkowania. Całkowanie sygnału w ruchomym oknie można uzyskać stosując wzór \eqref{eq:g}.
\begin{equation}
y(n) = \frac{1}{N}[x(n-(N-1)+)x(n-(N-2)+...+ x(n)]\label{eq:g}
\end{equation}
Na rysunku \eqref{qrsgoog} możemy zobaczyć rezultat prawidłowego doboru długości okna, cały zespół $QRS$ jest widoczny bez dodatkowych zbędnych informacji, na rysunku \eqref{qrsshor} widać efekt zbyt krótkiego okna całkującego, kompleks $QRS$ kończy się już przed załamkiem $R$, na rysunku \eqref{qrslong} widać efekt zbyt długiego okna całkowania, oprócz zespołu $QRS$ widać jeszcze część sygnału poza nim.
\begin{figure}
\centering
\includegraphics[scale=0.3]{include/qrsgood.eps}
\caption{Prawidłowo wykryty kompleks QRS (100\_ MLII.dat)}
\label{qrsgoog}
\end{figure}
\begin{figure}
\centering
\includegraphics[scale=0.3]{include/qrsshort.eps}
\caption{Za krótki wykryty kompleks QRS (100\_ MLII.dat)}
\label{qrsshor}
\end{figure}
\begin{figure}
\centering
\includegraphics[scale=0.3]{include/qrslong.eps}
\caption{Za długi wykryty kompleks QRS (100\_ MLII.dat)}
\label{qrslong}
\end{figure}
\newpage
\indent Postać sygnału otrzymanego po całkowaniu ruchomym oknem jest przedstawiona na rysunku \eqref{movingwindow}.\\
\begin{figure}
\centering
\includegraphics[scale=0.5]{include/movingwindow.eps}
\caption{Sygnał EKG scałkowany ruchomym oknem (100\_ MLII.dat)}
\label{movingwindow}
\end{figure}\\
\paragraph{Progowanie i detekcja zespołów QRS}
Ostatnim krokiem metody |verb|PanTompkins jest detekcja progu detekcji oraz sama detekcja zespołu QRS. Próg detekcji może być ustalony na stałe bądź wyliczony według własnych wzorów.
\paragraph{Wykrywanie załamka R}
\indent Aby odnaleźć załamek R należy teraz odnaleźć początek i koniec każdego kompleksu QRS w sygnale oryginalnym, używamy do tego punktów zapisanych w odnalezionym wektorze po operacji progowania. Mając taki wycinek oryginalnego sygnału, wyszukujemy wartości maksymalnej sygnału i oznaczamy ten punkt jako załamek R.\\
Na rysunku \eqref{rdet} można zobaczyć przykładową detekcję załamków R.
\begin{figure}[ht]
\centering
\includegraphics[scale=0.5]{include/rpeaks.eps}
\caption{Detekcja załamków R (100\_ MLII.dat)}
\label{rdet}
\end{figure}\\

\paragraph{Hilbert}
\label{sec:rs:papers}
\subparagraph{Usuwanie stałej DC i normalizacja sygnału}
Procedura operacji usuwania stałej DC i normalizacji sygnału przebiega w identyczny sposób jak opisano w procedurze "Pan-Thompkins"
\subparagraph{Filtr dolnoprzepustowy}
Filtracja dolnoprzepustowa przebiega w identyczny sposób jak opisano w procedurze "Pan - Thompkins"
\subparagraph{Filtr górnoprzepustowy}
Filtracja górnoprzepustowa przebiega w identyczny sposób jak opisano w procedurze "Pan - Thompkins"
\subparagraph{Różniczkowanie sygnału}
Po filtracji sygnału EKG, następnym krokiem jest różniczkowanie sygnału. Efektem takiego różniczkowania jest podkreślenie w sygnale nachyleń zespołów $QRS$ oraz stłumienie załamków $P$ i $T$. Pierwszą pochodną dla danego sygnału EKG można uzyskać za pomocą równania:
\begin{equation}
y(n) = \frac{1}{2 \Delta t}[x(n+1)-x(n-1)], dla n=0,1,2,3...,m-1\label{eq:h}
\end{equation}
gdzie m jest całkowita liczbą próbek w sygnale, a $\Delta t$ jest częstotliwością próbkowania.
Początkowy stan jest określony przez $x(-1)$ kiedy n = 0, a końcowy $x(m)$, kiedy n=m-1.
Warunki te minimalizują błąd na granicach.

\subparagraph{Podział sygnału w przesuwnym oknie}
Opisany poniżej algorytm transformaty Hilberta działa dobrze z krótkimi ciągami. Dlatego używa się okna o szerokości 1024 próbek do podziału zróżniczkowanego sygnału $y(n)$ przed uzyskaniem jego transformaty Hilberta. DO optymalizacji dokładności, początkowy punkt kolejnego okna powinien być zlokalizowany w punkcie ostatniego załamka $R$ w~poprzednim oknie. 
\subparagraph{Transformata Hilberta}
Transformata Hilberta rzeczywistej funkcji x(t) jest zdefiniowana jako:
\begin{equation}
\hat{x}=H[x(t)]=\frac{1}{\pi}\int_{-\infty}^\infty x(\tau)\frac{1}{t -\tau}d\tau \label{eq:i}
\end{equation}
Transformata Hilberta może być zinterpretowana jako konwolucja pomiędzy $x(t)$ i $\frac{1}{\pi t}$\\
Stosując transformatę Fouriera do \eqref{eq:i}, otrzymujemy:
\begin{equation}
F\lbrace \hat{x}(t)\rbrace= \frac{1}{\pi} F \lbrace \frac{1}{t} \rbrace F \lbrace x(t) \rbrace \label{eq:j}
\end{equation}
Ponieważ 
\begin{equation}
F\lbrace \frac{1}{t} \rbrace = \int_{-\infty}^{+\infty} \frac{1}{x} e^{-j 2 \pi fxdx} = -j \pi sgn f  \label{eq:k}
\end{equation}
gdzie $sgn f$ równe jest $+1$ dla $f > 0$, $0$ dla $f = 0$ i $-1$ dla $f<0$,\\
wtedy transformacja Fouriera transformaty Hilberta $x(t)$ jest podana przez \eqref{eq:j} jako:
\begin{equation}
F \lbrace \hat{x}(t) \rbrace= -j sgn f F \lbrace x(t) \rbrace \label{eq:l}
\end{equation}
W dziedzinie częstotliwości, wynik jest wtedy uzyskiwany przez przemnożenie widma $x(t)$ przez $j(+90^o)$ dla ujemnych częstotliwości i $-j(-90^o)$ dla dodatnich częstotliwości. Dziedzina czasowa może zostać uzyskana przez wyliczenie odwrotnej transformaty Fouriera. Dlatego, transformata Hilberta oryginalnej funkcji $x(t)$ reprezentuje jesgo harmoniczną koniugację.\\
Transformatę Hilberta $h(n)$ ciągu $y(n)$, który reprezentuje pierwszą pochodną sygnału EKG, jest uzyskiwana w według następujących kroków :
\begin{itemize}
		\item Uzyskaj transformatę Fouriera $F(n)$ wejściowego sygnału $y(n)$,
		\item Przemnóż dodatnie oraz ujemne harmoniczne odpowiednio przez $-j$ i $j$,
		\item Wyzeruj stałą składową $DC$,
		\item Wylicz odwrotną transformatę Fouriera poprzedniego ciągu.
\end{itemize}
\subparagraph{Adaptacyjny próg}
Szczyty w przetransformowanym sygnał metodą Hilberta $h(n)$ reprezentują obszary o wysokim prawdopodobieństwie wystąpienia prawdziwych załamków QRS. W praktyce, te szczyty różnią się od prawdziwych szczytów $R$ o kilka milisekund. W celu zagwarantowania dokładności detekcji załamków $R$, detektor drugiego poziomu jest wymagany. Ponieważ fale $P$ i $T$ są zminimalizowane w relacji do odpowiednich załamków w zespole QRS w ciągu Hilberta, detekcja prostym progiem jest używana do zlokalizowania załamków w ciągu $h(n)$. Próg musi być adaptacyjny w~celu zagwarantowania dokładności detekcji załamków R.\\
poziom progu jest ustalony według poniższych kryteriów:
\begin{itemize}
		\item Początkowo poziom szumu w analizowanym podciągu jest określany równoważną wartością $RMS(Root Mean Square)$ ciągu $h(n)$ i jego maksymalnej amplitudy w oknie, 
		\item Jeżeli wartość $RMS$ jest równa lub większa $18\%$ maksymalnej wartości ciągu $h(n)$, poziom szumu w segmencie uważany jest za wysoki i dlatego poziom progu ustawiany jest na $39\%$ maksymalnej amplitudy ciągu $h(n)$,  
		\item Jeżeli maksymalna wartość rozważanego ciągu jest dwukrotnie większa od maksymalnej wartości amplitudy poprzedniego okna 1024 próbek, wtedy próg jest powiększany do $39\%$ maksymalnej amplitudy poprzedniego okna $h(n)$, 
		\item Kiedy wartość $RMS$ jest niższa niż $18\%$ maksymalnej wartości ciągu $h(n)$, poziom szumu uważany jest za niski i poziom progu zmniejszany o $1.6$ raza wartości $RMS$,
		\item Jeżeli dwa załamki $R$ w ciągu Hilberta $h(n)$ są zlokalizowane blisko siebie(w mniejszym odstępie czasowym niż $200ms$), tylko jeden załamek jest prawdopodobnie rzeczywistym załamkiem $R$,
		\item Decyzja jest podejmowana na podstawie amplitudy załamka i ich pozycji w stosunku do ostatniego załamka $R$ zlokalizowanego za pomocą adaptacyjnego progu czasowego średniej długości interwału $R-R$ od poprzedniego zlokalizowanego załamka $R$.
\end{itemize}

\subsubsection{Opis procedur i metod}
\label{sec:rs:procs}

Implementacja modułu wykrywania załamków \verb|R| znajduje się w klasie \verb|RPeaksDetector| która rozszerza abstrakcyjny moduł wykrywania załamków \verb|R|. Klasa ta posiada implementację dwóch metod wykrywania załamków \verb|R|:
\begin{itemize}
	\item Pan-Tompkins
	\item Hilbert
\end{itemize}
Moduł na podstawie parametrów jakie otrzymuje z \verb|GUI| wybiera odpowiednią metodę detekcji. Oprócz metody detekcji możliwe jest również ustawienie ręczne parametrów takich jak:
\begin{itemize}
	\item dla metody Pan-Tompkins
	\begin{itemize}
		\item próg detekcji
		\item szerokość okna całkowania
	\end{itemize}
	\item dla metody Hilbert
	\begin{itemize}
		\item TODO JAKIE PARAMETRY
	\end{itemize}
\end{itemize}
Użytkownik może też skorzystać z automatycznej detekcji powyższych parametrów. Automatyczna detekcja jest włączona jako standardowe parametry.
Opis implementacji metod wykrywania załamków:
\paragraph{PanTompkins}
\begin{lstlisting}
bool panTompkinsRPeaksDetection(ECGSignalChannel *signal);
\end{lstlisting}
Funkcja otrzymuje na wejściu przefiltrowany sygnał \verb|ECGSignalChannel *signal| z modułu \verb|ECGBaseLine|. Z założenia sygnał jest pozbawiony składowej stałej, szumów pochodzących z mięśni orz znormalizowany w zakresie od -1 do 1. Detekcja pomija więc wstępną filtracje która została wykonana w poprzednim module. Zadania wykonywane wewnątrz można podzielić na następujące etapy:
\begin{itemize}
	\item różniczkowanie sygnału
	\item potęgowanie sygnału
	\item obliczanie szerokości okna całkowania 
	\item całkowanie sygnału ruchomym oknem
	\item obliczanie progu detekcji 
	\item wykrywanie zespołów QRS
	\item wykrywanie załamków R
\end{itemize}
Jeśli podczas wykrywania nie wystąpi żaden problem zostanie stworzony wektor z punktami w jakich zostały wykryte załamki \verb|R|. Natomiast jeśli z jakiegoś powodu wykrycie załamków nie powiedzie się zostanie rzucony wyjątek \verb|RPeaksDetectionExveption|.
\paragraph{Hilbert}
\begin{lstlisting}
bool hilbertRPeaksDetection(ECGSignalChannel *signal);
\end{lstlisting}
Funkcja również na wejściu otrzymuje sygnał z modułu \verb|ECGBaseLine|. Zadania wykonywane podczas detekcji można podzielić następująco:
\begin{itemize}
	\item różniczkowanie sygnału
	\item podział sygnału ruchomym oknem
	\item transformacja Hilberta
	\item wyznaczanie adaptacyjnego progu
	\item wykrywanie załamków R
\end{itemize}
Jeśli podczas wykrywania nie wystąpi żaden problem zostanie stworzony wektor z punktami w jakich zostały wykryte załamki \verb|R|. Natomiast jeśli z jakiegoś powodu wykrycie załamków nie powiedzie się zostanie rzucony wyjątek \verb|RPeaksDetectionExveption|

\subsubsection{Warunki testowania}
\label{sec:rs:tests}
Początkowo moduł był testowany przy użyciu własnego przefiltrowanego sygnału testowego utworzonego w Matlabie. Testowy sygnał został pozbawiony składowej stałej, zostały usunięte zakłócenia filtrami: dolnoprzepustowym oraz górnoprzepustowym. Cały sygnał został znormalizowany do zakresu od -1 do 1.\\
\indent Sygnał można wczytać przy pomocy funkcji jaka znajduje się w klasie \verb|RPeaksDetector| o nazwie \verb|getMockedSignal|. W pliku \verb|RPeaksDetector.h| znajdują się makrodefinicje które w prosty sposób pozwalają włączyć debugowanie procesu wykrywania sygnału oraz użycie testowego sygnału.\\

\indent Makrodefinicja \verb|USE_MOCKED_SIGNAL| włącza użycie testowego sygnału. Sygnał jaki dostarcza moduł \verb|ECGBaseLine| jest ignorowany. Dodatkowo makrodefinicja \verb|DEBUG| pozwala na wypisanie podstawowych informacji z przebiegu wykrywania załamków \verb|R|. Poniżej przedstawiamy przykładowy log zarówno dla metody \verb|PanTompkins| jak i \verb|Hilbert|:
\begin{verse}
Thersold size not found, use automatic calculated value\\ 
Input parameters for R peaks module:\\
Detection method: PanTompkins\\ 
Moving window size: 0\\ 
Thersold size: 0 \\
R peaks module started \\
Use mocked signal for R-peaks module. \\
Running module with custom parameters \\
Convolution [-0.125 -0.25 0.25 0.125] \\
Orginal signal size: 600000 \\
Exponentiation \^2 \\
Signal size after convolution: 599996 \\
Moving window integration\\
Calculating moving window size\\
Moving window size: 24 \\
Signal size after exponentiation: 599996 \\
Calculating detection thersold \\
After moving window integration signal size: 599972 \\
Final max value: 0.0438604 \\
Final mean value: 0.00221016\\
Current thresold value: 0.0109651 \\
Looking for points over thersold \\
Detect begin and end of QRS complex \\
Number of left points: 2093 \\
Number of right points: 2093 \\
Final R peaks detection \\
Number of detected R-peaks: 2093 \\
Done
\end{verse}

\begin{verse}
R peaks module started \\
Running module with default parameters \\
Window size not found, use automatic calculated value \\
Thersold size not found, use automatic calculated value \\
Input parameters for R peaks module:\\
Detection method: Hilbert \\
R peaks module started \\
Running module with custom parameters \\
Number of detected R-peaks: 2272 
\end{verse}


\indent Kiedy moduł \verb|ECGBaseLine| dostarczył nam prawidłowo przygotowany sygnał, rozpoczęliśmy testowanie na plikach pochodzących z bazy \verb|MIT_BIH|. Wyniki jakie otrzymaliśmy można znaleźć w kolejnym rozdziale.

\subsubsection{Wyniki}
\label{sec:rs:results}
Uzyskane wyniki wykrywania załamków \verb|R| dla sygnałów pochodzących z bazy \verb|MIT_BIH|.\\
\linebreak
\begin{tabular}{|c|c|c|c|c|c|}
\hline
sygnał & ilość R & \parbox{2.7 cm}{\centering wykryte \\ PanTompkins} & wykryte  Hilbert & \parbox{3 cm}{\centering skuteczność \\ PanTompkins} & skuteczność Hilbert\\ \hline
100   & 2273 & 3001 & 2272  & 91,7\% & 99,9\% \\ \hline
105   & 2572 & 2982 & 2559  & 85,1\% & 99,5\% \\ \hline
107   & 2137 & 3803 & 2125  & 33,1\% & 99,4\% \\ \hline
109   & 2532 & 2904 & 2492  & 85,4\% & 98,4\% \\ \hline
111   & 1774 & 3535 & 2089  & 12,2\% & 83,3\% \\ \hline
200   & 2601 & 3716 & 2593  & 66,8\% & 99,7\% \\ \hline
202   & 2136 & 3711 & 2116  & 36,3\% & 99\%   \\ \hline
220   & 2955 & 2055 & 2047  & 70\%   & 69,2\% \\ \hline
223   & 2332 & 3187 & 2501  & 63,2\% & 93,2\% \\ \hline
test100   & 2273 & 2274 & 2272  & 99,9\% & 99,9\% \\ \hline
\end{tabular}
\linebreak
\\
\indent Uzyskane wyniki pomiarów wyraźnie pokazują przewagę metody \verb|Hilberta|. W większości plików wynik skuteczności wahał się w około 98\%. Metoda detekcji \verb|PanTompkins| sprawia problemy głównie w~sygnałach o dużych zakłóceniach oraz tam gdzie morfologia pochodzenia sygnału uderzenia jest inna niż zatokowo-przedsionkowa.\\
\indent Ostatni sygnał 'test100' został wygenerowany z Matlaba, wczytaliśmy sygnał 100.dat, przefiltrowaliśmy, usunęliśmy składową stałą oraz znormalizowaliśmy. Tak przygotowany sygnał został wczytany i przeprowadziliśmy detekcję. Jak widać skuteczność jest znacznie większa niż w przypadku oryginalnego sygnału gdzie przygotowany sygnał dostajemy z modułu \verb|BaseLIne|.\\
\indent Próbowaliśmy również testować wykrywanie zmieniając ostawienia modułu \verb|BaseLine| wybierając inne filtracje. Wyniki jakie uzyskaliśmy były praktycznie identyczne jak powyżej dlatego nie zamieszczamy ich tutaj.

\subsection{Waves}
\label{sec:waves}
Autorzy: Agata Sitnik i Łukasz Zieńkowski.

\subsubsection{Opis zadania}
\label{sec:waves:desc}
Celem modułu było wykrycie zespołów QRS oraz załamków P i T.
Zespół QRS to największy zespół załamków EKG. Opisuje depolaryzację (pobudzenie) mięśni komór serca i składa się z jednego lub kilku załamków określanych kolejno jako Q, R i S, w zależności od miejsca wystąpienia i kierunku wychylenia. Czas trwania zespołu prawidłowo wynosi od 0,6-0,11 s i składa się z następujących składowych:
\begin{itemize}
     \item \verb|Załamek Q|  -- pierwsze ujemne wychylenie zespołu QRS
     \item \verb|Załamek R|  -- pierwsze dodatnie wychylenie zespołu QRS
     \item \verb|Załamek S|  -- każde ujemne wychylenie za załamkiem R
\end{itemize}

Poza zespołem QRS możemy wyróżnić również załamki P i T.
Załamek P odpowiada depolaryzacji przedsionków. W warunkach prawidłowych jest on dodatni w odprowadzeniach I, II, aVF i ujemny w~aVR; czas jego trwania jest krótszy od 0.12 s, a amplituda nie przekracza 2.5 mm w odprowadzeniach kończynowych.
Odcinek PQ odpowiada repolaryzacji przedsionków. W warunkach prawidłowych repolaryzacja przedsionków nie powoduje przemieszczenia odcinka PQ w stosunku do odcinka TP i~odcinek PQ przebiega w linii izoelektrycznej.
Odstęp PQ jest elektrokardiograficznym odpowiednikiem czasu wędrowania bodźca z węzła zatokowego
przez prawy przedsionek, węzeł przedsionkowo - komorowy, pęczek Hisa i włókna Purkinjego aż do mięśnia komór. W warunkach prawidłowych czas jego trwania zależy od wieku badanej osoby i częstotliwości rytmu serca mieszcząc się w granicach od 0.12 s do 0.20 s; u osób w starszym wieku za górną granicę normy można przyjąć wartość 0.23 s.
Odcinek ST jest elektrokardiograficznym odpowiednikiem czasu wędrowania bodźca z węzła zatokowego przez prawy przedsionek, węzeł przedsionkowo - komorowy, pęczek Hisa i włókna Purkinjego aż do mięśnia komór. W warunkach prawidłowych czas jego trwania zależy od wieku badanej osoby i częstotliwości rytmu serca mieszcząc się w granicach od 0.12 s do 0.20 s; u osób w starszym wieku za górną granicę normy można przyjąć wartość 0.23 s.
Załamek T odpowiada fazie szybkiej repolaryzacji mięśnia komór. Prawidłowy załamek T jest dodatni w odprowadzeniach I, II, aVL, V3 - V6 a ujemny w aVR.


\subsubsection{Badania literaturowe}
\label{sec:waves:papers}
W literaturze można znaleźć kilka rozwiązań problemu odnalezienia zespołów QRS i załamków P i T. Większość z nich opiera się na dyskretnej transformacie falkowej, która polega na splocie
sygnału pochodzącego z zapisu badania EKG z wybraną falką, pobraną z bazy. Baza jest stale rozwijana o nowe falki. Filtry cyfrowe analizy falkowej odpowiadają falce i funkcji skalującej w postaci tzw.
spline’ów bramkowych drugiego stopnia o zwartym i krótkim nośniku. Dzięki temu
podczas analizy sygnału i~detekcji osobliwości możemy dokładniej kontrolować parametry
procesu separacji wybranych częstotliwości. Dzięki analizie wieloskalowej możliwe jest
zlokalizowanie miejsca gwałtownej zmiany sygnału, a tym samym lokalizacji zespołu QRS.
Metoda posiada mniejszą wrażliwość na zmiany morfologii kolejnych zespołów QRS,
minimalizuje problemy związane z występowaniem składowej wolnozmiennej, artefaktów
ruchu i napięcia mięśni oraz pozwala na łatwiejszą separację załamka R w stosunku do
załamków P i T.

Wiele rozwiązań zagadnienia detekcji zespołów QRS opiera się na filtracji adaptacyjnej,
zastosowaniu ukrytych modeli Markowa, algorytmów genetycznych lub
transformacie Hilberta. Wszystkie opierają się na pomyśle utworzenia detektora zespołów QRS, jednak pozwalają jedynie na drobną poprawę
skuteczności lub szybkości działania, gdyż tak naprawdę nie znaleziono jeszcze w 100\% skutecznej funkcji detekcyjnej. 
Istnieją dane literaturowe, w których detekcja zespołów QRS opiera się na stosowaniu filtrów wysokopasmowych. Filtry pełnią rolę narzędzia różniczkującego, a wykorzystywane są pierwsze i drugie pochodne sygnału. Różniczkowanie w tych metodach wykonuje się w celu uwypuklenia cech charakterystycznych sygnału. Zróżniczkowany sygnał poddaje się progowaniu. Pozwala to na wykrycie potencjalnych zespołów QRS, a dodatkowe progi pozwalają eliminować fałszywych kandydatów. Próg składa się z trzech parametrów: adaptacyjnej wartości slew-rate, drugiej wartości, która rośnie gdy w sygnale obserwuje się wzrost częstotliwości i trzeciej mającej na celu uniknięcie braku niskich wartości amplitudy. Zaletą tych metod jest ich prostota, a wadą niska efektywność i skomplikowane ciągi decyzyjne. Istnieją również metody łączące tosowanie filtrów wysokoprzepustowych z~niskoprzepustowymi, pozwala to zwiększyć efektywność wykrywania zespołów QRS.


Jedną z najbardziej innowacyjnych metod automatycznej detekcji zespołów QRS
jest stosowanie do tego celu sieci neuronowych typu wielowarstwowego
perceptronu lub sieci LVQ. Budowa sieci MLP opiera się na konstrukcji z prostych neuronów o funkcji przejścia w postaci
tangensa hiperbolicznego lub sigmoidy. Liczba warstw sieci nie powinna
przekraczać 3, a liczba neuronów również powinna być ograniczona. 
Sieć neuronowa MPL działa na zasadzie próby przewidzenia stanu
przyszłego w oparciu o stan poprzedni. Wadą tego rozwiązania jest to, że najpierw należy przeprowadzić uczenie sieci oraz duża wrażliwość na zaszumienie sygnału.  Sieci typu LVQ działają inaczej: neurony środkowej warstwy uczą się rozpoznawania
sygnałów określonego typu. Ich liczba określa liczbę grup sygnałów.
Wadą sieci LVQ jest to, że jest bardzo rozbudowana oraz
wymaga dosyć dużego zbioru testowego. Jednak dobrze nauczona sieć typu LVQ działa w zasadzie dla każdego zapisu EKG.
Obie metody MLP i LVQ wymagają jednak
dużego nakładu pracy przy tworzeniu sieci oraz ich uczeniu. Ważne jest również to, ze błędnie przeprowadzony sposób uczenia, może znacząco pogorszyć wykrywalność zespołów QRS.

W przypadku tego projektu zdecydowano się na wykorzystanie algorytmu treshold based detection.


\subsubsection{Opis procedur i metod}
\label{sec:waves:procs}
Detekcja opiera się na algorytmie \verb|treshold based detection|.


zastosowane funkcje:


\begin{lstlisting}
bool QRSPointsDetector::detectQRS()
\end{lstlisting}
Główna funkcja odpowiedzialna za detekcję zespołów QRS. Funkcja korzysta ze znajomości położenia
załamków R, uzyskanych z modułu Rpeaks. Określa ona minima po obu stronach tych załamków,
następnie wykonuje normalizację sygnału, na którą składa się:
\begin{itemize}
     \item odjęcie średniej z całego sygnału
     \item spotęgowanie (wzmocnienie)
     \item podzielenie przez maksymalną wartość występującą w sygnale (otrzymujemy sygnał
w zakresie 0-1)
\end{itemize}

Po normalizacji w funkcji jest wydzielana część sygnału o intensywności powyżej 5\%, które ma na
celu wyzerowanie sygnału poza obszarem QRS. Funkcja służy do obliczenia parametrów \verb|QRS_onset| i~\verb|QRS_end|.

\begin{lstlisting}
bool QRSPointsDetector::detectPT()
\end{lstlisting}

Główna funkcja odpowiedzialna za detekcję załamków P i T. Funkcja najpierw wstępnie przetwarza
otrzymany sygnał, wraz z oznaczonymi miejscami występowania zespołów QRS. Pierwszym etapem jest
wyzerowanie sygnału w miejscach występowania QRS. Następnie wykonywana jest normalizacja sygnału
składająca się z następujących etapów:
\begin{itemize}
     \item Każdą część sygnału znajdującą się w okresie podzielić przez maksymalną wartość w danym okresie
     \item Filtracja sygnału
     \item Odjęcie mediany dla danego okresu od wszystkich wartości w okresie
     \item wydzielenie części sygnału o intensywności powyżej 7\% (usunięcie sygnału poza
QRS)
\end{itemize}

Funkcja następnie wykonuje właściwą detekcję załamków P i T. Ich wykrycie opiera się na następujących
założeniach:
\begin{itemize}
     \item Załamek T zawsze jest dłuższy niż 90 próbek (najmniejszy zanotowany to 95, a
największy 208).
     \item Długość załamka jest większa niż 30\% długości okresu.
     \item Odległość między końcem załamka T, a początkiem następnego QRS jest zawsze
większa niż 25\% długości okresu.
     \item Załamek P zawsze jest dłuższy niż 9 próbek (najmniejszy zanotowany to 10).
     \item Odległość między początkiem załamka P, a początkiem następnego QRS jest
zawsze mniejsza niż 23\% długości okresu.
\end{itemize}

W celu wyznaczenia załamków P i T w funkcji są najpierw wyznaczane początki załamków T jako
końce poprzedzającego go zespołu QRS . W danym okresie sygnału wynikowego, funkcja znajduje
pierwszą próbkę, której wartość jest różna od zera (dalej pod nazwą próbka pierwsza). Następnie
znajduje pierwszą próbkę po tej poprzedniej, której wartość jest równa zero (dalej pod nazwą próbka
druga). Jeżeli odległość pomiędzy początkiem następnego zespołu QRS, a tą pierwszą próbką
jest większa niż 25\% długości okresu i jednocześnie odległość między drugą próbką, a końcem
poprzedzającego ją QRS jest większa niż 90, oznacza to, że druga próbka jest końcem załamka T.
Pierwsze z tych założeń zapobiega zidentyfikowaniu załamka P jako T, w przypadku gdy T jest
odwrócone i~w~sygnale wynikowym nie występuje. Drugie z założeń mówiące o szerokości załamka
T, zapobiega zaklasyfikowaniu szumów o małej długości jako załamka T.


\begin{lstlisting}
double findMinimum (ECGSignalChannel *signal,int forBegin, int forEnd)
\end{lstlisting}
Funkcja znajdująca minimum lokalne w danym zakresie sygnału, parametry forBegin i forEnd
określają początek i koniec analizowanego okna.

\begin{lstlisting}
double findMaximum (ECGSignalChannel *signal,int forBegin, int forEnd)
\end{lstlisting}
Funkcja znajdująca maksimum lokalne w danym zakresie sygnału, parametry forBegin i forEnd
określają początek i koniec analizowanego okna.

\begin{lstlisting}
void QRSPointsDetector::runModule(const ECGSignalChannel &filteredSignal,
 const ECGInfo &ecgi, const ECGRs &ecgRs, ECGWaves & ecgWaves)
\end{lstlisting}
Funkcja uruchamiająca moduł

Parametry:
\verb|&filteredSignal| -- przefiltrowany sygnał

\verb|& ecgi| -- informacje o sygnale

\verb|&ecgRs| -- położenie załamków R

\verb|& ecgWaves| -- informacja zwrotna -- położenie zespołów QRS oraz załamków P i T

\begin{lstlisting}
ECGSignalChannel QRSPointsDetector::gradient(ECGSignalChannel * signal)
\end{lstlisting}

\begin{lstlisting}
ECGSignalChannel QRSPointsDetector::averageFilter(ECGSignalChannel * signal)
\end{lstlisting}

Funkcja realizująca filtr uśredniający sygnał

\begin{lstlisting}
ECGRs QRSPointsDetector::getMockedRPeak()
\end{lstlisting}
Funkcja ustawiająca wartości załamków R, używana do testowania modułu.

\begin{lstlisting}
ECGSignalChannel QRSPointsDetector::getMockedSignal()
\end{lstlisting}
Funkcja ustawiająca wartości sygnału, używana do testowania modułu.



\subsubsection{Warunki testowania}
\label{sec:waves:tests}

W celu przetestowania algorytmu najpierw wykonano implementację w programie Matlab 2009, sprawdzenie skuteczności algorytmu wykonano na sygnałach ściągniętych z bazy wykorzystywanej w~poprzednim semestrze na zajęciach z przedmiotu Przetwarzanie sygnałów w systemach diagnostyki medycznej. 
Następnie kod przepisano w C++ w środowisku VisualStudio 2010 i tam testowano go za pomocą tych samych sygnałów. Posłużono się jako danymi wejściowymi również wygenerowanym wektorem położenia załamków R.

\subsubsection{Wyniki}
\label{sec:waves:results}
\begin{center}
%
\begin{figure}
\begin{centering}
\includegraphics[scale=0.4]{include/waves_matlab}
\par\end{centering}

\caption{Wyniki modułu Waves -- Matlab, oznaczone zespoły QRS i załamki R}
\label{fig:waves_matlab}
\end{figure}

\par\end{center}


\begin{center}
%
\begin{figure}
\begin{centering}
\includegraphics[scale=0.4]{include/waves_matlab2}
\par\end{centering}

\caption{Wyniki modułu Waves -- Matlab, oznaczone załamki P, T i R }
\label{fig:waves_matlab2}
\end{figure}

\par\end{center}

\begin{center}
%
\begin{figure}
\begin{centering}
\includegraphics[scale=0.4]{include/waves_cpp}
\par\end{centering}

\caption{Wyniki modułu Waves -- Visual Studio 2010}
\label{fig:waves_cpp}
\end{figure}

\par\end{center}

Udało się osiągnąć zadowalające wyniki, które nieznacznie odbiegają od danych przewidywanych. 
Wyniki osiągnięte podczas wstępnej implementacji w Matlabie, gdzie wykonano oznaczenia początków i końców zespołów QRS \ref{fig:waves_matlab} oraz następnie oznaczono położenie załamków P i końców załamków T \ref{fig:waves_matlab2} są nieznacznie gorszej jakości niż uzyskane później w wynikowym programie \ref{fig:waves_cpp} . Jest to spowodowane między innymi tym, że program w Matlabie był przygotowany jako prototyp, który później został dopracowany i rozszerzony o dodatkowe warunki klasyfikacji sygnału jako załamki oraz o dodatkową filtrację.

\subsection{HRV1}
\label{sec:hrv1}

Autorzy: Łukasz Jaromir i Leszek Sosnowski
\subsubsection{Opis zadania}
\label{sec:hrv1:desc}
\begin{description}
\item[Temat] Analiza zmienności rytmu serca (HRV) \#1
\item[Opis] Zmienność rytmu zatokowego opisuje różnice w długościach interwałów RR
wyznaczanych przez kolejne szczyty zespołów QRS. Występowanie tych różnic świadczy
o zdolności serca do adaptacji względem zewnętrznych bodźców i przedstawia informacje
o pracy autonomicznego systemu nerwowego. Najprostsze ze stosowanych metod to analiza
statystyczna w dziedzinie czasu oraz analiza częstotliwościowa. Celem projektu jest
opracowanie i implementacja metod czasowych oraz częstotliwościowych analizy zmienności
rytmu serca.
\item[Dane] ciąg próbek załamków R z modułu RPEAKS
\item[Szukane] moduł programu wyznaczający, wyświetlający parametry analizy czasowej
i częstotliwościowej, a także rysujący postać częstotliwościową tachogramu wraz
z naniesionymi zakresami parametrów HF, LF, VLF, ULF
\end{description}

\subsubsection{Badania literaturowe}
\label{sec:hrv1:papers}
Ocena zmienności rytmu serca jest podstawową techniką umozliwiającą ocenę wpływu układu współczulnego (pobudzanie) i przywspółczulnego (hamowanie) na akcję serca. Ponieważ rola tych układów silnie zależy od stanu człowieka (sen, czuwanie) dopiero analiza całodobowych zmian rytmu serca pozwala diagnozować poprawność tego współdziałania.

Analiza zmienności rytmu serca opiera się na bieżących wartościach interwałów międzyuderzeniowych (odstępów RR) i może być prowadzona:
\begin{itemize}
\item w dziedzinie czasu (metodami statycznymi),
\item metodami geometrycznymi (aproksymacja histogramu trójkątem),
\item metodami częstotliwościowymi (z wykorzystaniem funkcji autokorelacji lub przekształcenia Fouriera).
\end{itemize}

Tachogram (po uzupełnieniu luk po zespołach QRS innych niż nadkomorowe) może być podstawą wyznaczania statycznych i geometrycznych współczynników zmienności akcji serca:
\newline
\begin{tabular}{|l|p{14cm}|}
\hline
Oznaczenie: & Znaczenie wskaźnika \\ \hline
Średni RR & Wartość średnia ze wszystkich odstępów RR rytmu zatokowego \\ \hline
SDNN & odchylenie standardowe interwałów RR (ang. standard deviation normal-to-normal) \\ \hline
rMSSD & pierwiastek kwadratowy ze średniej kwadratów różnic pomiędzy kolejnymi dwoma interwałami RR \\ \hline
NN50 & liczba interwałów RR, których różnica przekracza 50 ms \\ \hline
pNN50 & odsetek różnic pomiędzy interwałami RR, które przekraczają 50 ms [\%] \\ \hline
SDANN & odchylenie standardowe ze wszystkich średnich interwałów RR w 5 minutowych segmentach czasu całego zapisu [ms] \\ \hline
SDANNindex & średnia z odchyleń standardowych interwałów RR w 5 minutowych segmentach czasu całego zapisu [ms] \\ \hline
SDSD & odchylenie standardowe różnic pomiędzy dwoma sąsiadującymi interwałami RR [ms] \\
\hline
\end{tabular}
\linebreak

Analiza zmienności rytmu serca w dziedzinie częstotliwości wymaga rozwiązania problemu dotyczącego natury sygnału: tachogram jest dyskretną funkcją czasu próbkowaną niejednorodnie (interwał próbkowania nie jest stały). Konieczne jest więc odtworzenie ciągłej funkcji czasu na podstawie posiadanego tachogramu, a następnie ponowna jej dyskretyzacja w sposób jednorodny.
Algorytm interpolacji, który wybraliśmy to algorytm interpolacji funkcjami sklejanymi trzeciego stopnia.

Kolejnym krokiem jest stosowanie przekształcenia Fouriera w celu otrzymania widma sygnału. Wykres widma dzielony jest zwykle na pasma, w których określana jest całkowita moc wyrażana następnie przy użyciu współczynników:
\newline
\begin{tabular}{|c|p{14cm}|}
\hline
Oznaczenie: & Znaczenie wskaźnika \\ \hline
TP & całkowita moc widma (większe lub równe 0.4 Hz)  \\ \hline
HF & moc widma w zakresie wysokich częstotliwości (0.15 - 0.4 Hz) \\ \hline
LF & moc widma w zakresie niskich częstotliwości (0.04 - 0.15 Hz) \\ \hline
VLF & moc widma w zakresie bardzo niskich częstotliwości (0.003 - 0.04 Hz)  \\ \hline
ULF & moc widma w zakresie ultra niskich częstotliwości (0.003 Hz)  \\ \hline
LFHF & stosunek mocy widm w zakresie niskich częstotliwości do wysokich częstotliwości  \\
\hline
\end{tabular}


Analiza częstotliwościowa dostarcza informacji o dynamicznej równowadze układu autonomicznego umożliwiając ocenę współpracy układu współczulnego i przywspółczulnego. Współczynniki analizy częstotliwościowej HRV dobrze reprezentują zjawiska zmienności \cite{AUGUST1}.

\subsubsection{Opis procedur i metod}
\label{sec:hrv1:procs}
Klasą realizującą obliczenia parametrów analizy czasowej i częstotliwościowej jest klasa HRV1Analyzer.
Metody w klasie HRV1Analyzer:
\begin{itemize}
\item \verb+HRV1Analyzer::HRV1Analyzer()+ -- konstruktor klasy
\item \verb+HRV1Analyzer::~HRV1Analyzer()+ -- destruktor klasy
\item \verb+void HRV1Analyzer::runModule(const ECGRs & r_peaks_data, ECGHRV1 & hrv1_data)+ -- metoda wywoływana przez kontroler służąca do uruchomienia modułu i realizacji wyliczenia współczynników
\item \verb+void HRV1Analyzer::prepareSignal()+ -- funkcja zajmująca się przygotowaniem sygnału -- metoda przekształca tablicę z numerami próbek na tablicę zawierającą czas pomiędzy 2 załamkami
\item \verb+void HRV1Analyzer::prepareSigAbsolute()+ -- funkcja zajmująca się dalszym przygotowaniem sygnału -- tworzy sygnał zawierający bezwzględny czas wystąpienia kolejnych załamków RR
\item \verb+void HRV1Analyzer::calculateParameters()+ -- właściwa metoda odpowiadająca za wyliczenie wszystkich współczynników ilościowych oraz częstotliwościowych, które składają się na właściwą analizę, za którą odpowiada moduł HRV1
\item \verb+void HRV1Analyzer::setParams(ParametersTypes &parameterTypes)+ -- metoda odpowiedzialna za ustawianie parametrów (w naszym module nie wykorzystywana)
\item \verb+ouble* HRV1Analyzer::doFFT(double* sigAfterSpline)+ -- metoda przeprowadzająca analizę FFT zinterpolowanego sygnału z wykorzystaniem biblioteki KISS FFT
\item \verb+kiss_fft_cpx* HRV1Analyzer::copycpx(double *mat, int nframe)+ -- funkcja alokująca klasę z danymi potrzebnymi do realizacji transformaty Fouriera
\item \verb+double* HRV1Analyzer::cubicSpline(double* x, double* y, int nframe)+ -- funkcja odpowiedzialna za interpolowanie próbek (wykonywana jest interpolacja funkcjami sklejanymi trzeciego stopnia). Zachodzi interpolacja funkcji o wektorach x i y  -- długość wektorów x i y jest taka sama i wynosi nframe
\item \verb+double HRV1Analyzer::mean(double *tab, int start, int end)+ -- funkcja obliczająca średnią arytmetyczną elementów z tablicy tab
\item \verb+double HRV1Analyzer::std(double *tab, int start, int end)+ -- funkcja obliczająca odchylenie standardowe
\end{itemize}

Klasą przechowującą współczynniki analizy częstotliwościowej i czasowej jest klasa ECGHRV1.

Szukając biblioteki odpowiedzialnej za realizację transformaty Fouriera skupiliśmy się na znalezieniu narzędzia maksymalnie prostego i szybkiego aby kod pozostał zoptymalizowany. Zdecydowaliśmy się na darmową bibliotekę KISS (Keep it Simple, Stupid).

Tymi samymi przesłankami kierowaliśmy się szukając biblioteki odpowiedzialnej za interpolowanie sygnału. Nasz wybór padł na bibliotekę algib.

\subsubsection{Warunki testowania}
\label{sec:hrv1:tests}

Jako, że realizacja projektu odbywała się w grupach i tempo pracy było różne musieliśmy sobie poradzić i rozpocząć testowanie nie mając R\_PEAKów. Do tego celu stworzyliśmy prosty algorytm, który emulował pracę całego kontrolera. Program ten na wejście podawał zmienne, które w przyszłości otrzymalibyśmy z grupy R\_PEAKs  -- korzystaliśmy z materiałów, które otrzymaliśmy w poprzednim semestrze realizując w Matlabie laboratorium HRV1. Program zwracał natomiast policzone współczynniki analizy czasowej i częstotliwościowej.

W ten sposób, mogliśmy skutecznie pracować nad optymalizacją algorytmów porównując otrzymane wyniki z tymi, które pochodziły z prototypu napisanego w Matlabie, czekając na sprzężenie naszego modułu z modułem R\_PEAKS.

\subsubsection{Wyniki}
\label{sec:hrv1:results}
Pierwszy etap prac, który zakończyliśmy implementując metody odpowiedzialne za obliczanie poszczególnych współczynników, pozwolił nam na ich porównanie do tych otrzymywanych z prototypu napisanego w środowisku Matlab. Były one zbliżone, ale nie satysfakcjonujące: brak algorytmu aproksymacji w sposób istotny zwiększał błędy.

Dopiero dalsze poprawki i optymalizacja, a przede wszystkim implementacja algorytmu aproksymacji skutecznie zminimalizował błędy, tak że wyniki przez nas otrzymane były zadowalające.


\subsection{HRV2}
\label{sec:hrv2}

Autor: Krzysztof Farganus


\subsubsection{Opis zadania}
\label{sec:hrv2:desc}

Celem projektu jest opracowanie i implementacja metod geometrycznych
analizy HRV.

Dane przyjmowane przez moduł: 
\begin{itemize}
\item ciąg próbek załamków R z modułu R\_PEAKS.
\end{itemize}
Dane zwracane przez moduł:
\begin{itemize}
\item wykres histogramu 
\item wskaźnik TINN 
\item Indeks Trójkątny
\item wykres Poincare wraz z parametrami SD1 i SD2
\end{itemize}

\subsubsection{Badania literaturowe}
\label{sec:hrv2:papers}

Techniki geometryczne służą do przedstawienia długookresowej zmienności
rytmu serca \cite{hrv2-pl}. Są łatwe do uzyskania, ponieważ bazują na aproksymacji
histogramu trójkątem. Szerokość przedziałów klasowych histogramu ma
tutaj kluczowe znaczenie, gdyż jej wartość wpływa na rezultaty metod
geometrycznych. Głównie stosowany jest zakres wynoszący 7.8125 ms
(1/128 s). Powodem jest częstotliwość próbkowania sygnału o najczęściej
występującej wartości 128 Hz.

Cechy metod geometrycznych: 
\begin{itemize}
\item odporność na zakłócenia ze względu na zastosowanie technik aproksymacyjnych
\item eliminacja artefaktów zlokalizowanych poza trójkątem 
\item wyniki niezależne od jakości zapisu sygnału 
\item rezultaty zależne od czasu rejestracji 
\item wymagana duża liczba odstępów RR dla poprawnej analizy (minimalny
czas zapisu \textendash{} 20 minut) 
\end{itemize}
Do najbardziej rozpowszechnionych metod należą: 
\begin{itemize}
\item Wykres histogramu przedstawiający rozkład interwałów RR
\item Indeks trójkątny (HRV triangular index )  -- całkowita liczba wszystkich
odstępów RR podzielona przez liczbę odstępów RR o najczęściej spotykanym
czasie trwania. 
\item Trójkątna interpolacja odstępów RR (TINN) \textendash{} długość podstawy
trójkąta aproksymującego histogram kolejnych odstępów interwałów RR
rytmu zatokowego wyrażana w milisekundach. 
\item Wykres Poincare  -- graficzna reprezentacja korelacji pomiędzy kolejnymi
interwałami , gdzie każdy odstęp RR jest opisany funkcją RR+1. Do
analizy rozproszenia punktów na wykresie stosuje się dwa deskryptory:
SD1 oraz SD2, odpowiadające odchyleniom standardowym. Pierwszy charakteryzuje
rozkład punktów w poprzek linii identyczności, natomiast drugi wzdłuż
tej linii.
\end{itemize}
\textbf{Algorytm obliczania wskaźnika TINN \cite{hrv2-eng}}

Aby obliczyć parametr TINN, czyli wyznaczyć wartości punktów N i M
należy opracować funkcję multiliniową q(t) o postaci:
\begin{itemize}
\item $q(t) = 0$ dla \ensuremath{t \le N}
\item $q(t) = 0$ dla \ensuremath{t \ge M} 
\item $q(X) = Y$ w pozostałych przypadkach
\end{itemize}
a następnie znaleźć minimum z całki o wzorze:

\begin{equation}
\int_{0}^{+\infty}\left(\: D(t)\right)-q(t)\:)^{2}dt
\end{equation}

spośród wszystkich kombinacji punktów (N,M) , gdzie D(t) to wartość
histogramu. W układzie dyskretnym poprzedni wzór wygląda następująco:

\begin{equation}
\sum(\: D(t)-q(t)\:)^{2}\rightarrow \text{minimum}
\end{equation}


przy czym dla $t\:\epsilon(0,N)$ oraz $t\:\epsilon(M,\infty)$ ma
postać: 

\begin{equation}
D(t)^{2}
\end{equation}


natomiast dla $t\:\epsilon\left\langle 0,N\right\rangle$ wygląda
następująco: 
\begin{equation}
(D(t)-q(t))^{2}
\end{equation}

Algorytm obliczania parametrów SD1 i SD2

Parametry SD1 i SD2 są wyznaczane według poniższych wzorów:

\begin{align}
SD1 &= \sqrt{\frac{1}{2}\cdot SDSD^{2}} \\
SD2 &= \sqrt{2\cdot SDNN^{2}-\frac{1}{2}\cdot SDSD^{2}}
\end{align}


gdzie:
\begin{itemize}
\item SDNN to odchylenie standardowe interwałów RR: $SDNN=\sqrt{\frac{1}{N-1}\cdot\sum_{i=1}^{N}\left(\bar{RR}-RR_{i}\right)^{2}}$
\item SDSD to odchylenie standardowe różnic pomiędzy dwoma sąsiadującymi
interwałami: $SDSD=\sqrt{E\left\{ \triangle RR_{j}^{2}\right\} -E\left\{ \triangle RR_{j}\right\} ^{2}}$
\end{itemize}

\subsubsection{Opis procedur i metod }
\label{sec:hrv2:procs}

\paragraph{\texttt{GeometricAnalysis::runModule}} -- wirtualna funkcja uruchamiająca moduł HRV2.

\medskip{}


Argumenty funkcji:
\begin{itemize}
\item \verb+ECGInfo info+ -- dane o wczytanym sygnale
\item \verb+ECGRs ecgRs+ -- wektor numerów próbek zawierający załamki R
\item \verb+ECGHRV2 ecgHRV2+ -- instancja klasy przechowującej wyniki analizy zmienności
rytmu serca metodami geometrycznymi
\end{itemize}
\medskip{}


Funkcja zwraca:

Funkcja przekazuje wyniki metod geometrycznych analizy HRV do instancji
klasy \verb+ECGHRV2+.

\medskip{}


Używane funkcje:
\begin{itemize}
\item \verb+PrepareRRSignal+
\item \verb+MakeHistogramAndGeometricalParams+
\item \verb+MakePoincareAndSDParams+
\item \verb+SetHRV2Params+\medskip{}

\end{itemize}
Używane zmienne:
\begin{itemize}
\item \verb+ECGRs rpeaks+ -- atrybut klasy \verb+GeometricAnalysis+ przechowujący wektor
numerów próbek z~załamkami R
\item \verb+double SamplingInterval+ -- atrybut klasy \verb+GeometricAnalysis+ przechowujący
częstotliwość wczytanego sygnału
\end{itemize}


\medskip{}

\paragraph{\texttt{GeometricAnalysis::PrepareRRSignal}} -- funkcja przekształca wektor numerów próbek z załamkami R na wektor
interwałów RR w milisekundach.

\medskip{}


Argumenty funkcji:
\begin{itemize}
\item \verb+rpeaks+ -- wektor numerów próbek z załamkami R
\end{itemize}
\medskip{}


Funkcja zwraca:

Funkcja zapisuje wektor interwałów RR w atrybucie klasy \verb+RR_intervals+.

\medskip{}


Używane funkcje:
\begin{itemize}
\item \verb+gsl_vector_int_get+ -- metoda GSL pobierająca całkowitą wartość danego
elementu z wektora
\item \verb+gsl_vector_set+ -- metoda GSL zapisująca zmiennoprzecinkową wartość
w danym elemencie wektora
\item \verb+gsl_vector_scale+ -- metoda GSL mnożąca każdy element wektoru przez
liczbę zmiennoprzecinkową
\end{itemize}
\medskip{}


Używane zmienne:
\begin{itemize}
\item \verb+unsigned int rpeaks_size+ -- długość wektora z numerami próbek załamków
R
\item \verb+double SamplingInterval+ -- częstotliwość wczytanego sygnału
\end{itemize}

\medskip{}

\paragraph{\texttt{GeometricAnalysis::MakeHistogramAndGeometricalParams}} -- funkcja tworzy wektor wartości histogramu, oblicza wysokość i pozycję
kolumny histogramu reprezentującą najczęściej powtarzający się interwał
RR, wylicza długość oraz pozycję podstawy trójkąta aproksymującego,
oszacowuje wartość indeksu trójkątnego.

\medskip{}


Argumenty funkcji:
\begin{itemize}
\item \verb+RR_intervals+ -- wektor interwałów RR w milisekundach
\end{itemize}
\medskip{}


Funkcja zwraca:
\begin{itemize}
\item histogram\_x -- pozycje kolumn histogramu
\item histogram\_y -- wysokości kolumn histogramu
\item X -- pozycja najwyższej kolumny histogramu
\item Y -- wysokość najwyższej kolumny histogramu
\item N -- początek podstawy trójkąta aproksymującego
\item M -- koniec podstawy trójkąta aproksymującego
\item HRVTriangularIndex -- indeks trójkątny
\item TINN -- wskaźnik TINN
\end{itemize}
\medskip{}


Używane funkcje:
\begin{itemize}
\item \verb+gsl_vector_max+ -- metoda GSL zwracająca największy element danego
wektora
\item \verb+gsl_vector_min+ -- metoda GSL zwracająca najmniejszy element danego
wektora
\item \verb+gsl_vector_int_max_index+ -- metoda GSL zwracająca pozycję największego
elementu danego wektora
\item \verb+gsl_vector_int_get+ -- metoda GSL pobierająca całkowitą wartość danego
elementu z wektora
\item \verb+gsl_vector_get+ -- metoda GSL pobierająca zmiennoprzecinkową wartość
danego elementu z~wektora
\item \verb+gsl_vector_set+ -- metoda GSL zapisująca zmiennoprzecinkową wartość
w danym elemencie wektora
\item \verb+gsl_vector_int_set+ -- metoda GSL zapisująca całkowitą wartość w
danym elemencie wektora
\item \verb+gsl_interp_alloc+ -- metoda GSL tworząca wskaźnik do nowo utworzonego
obiektu interpolacji
\item \verb+gsl_interp_accel_alloc+ -- metoda GSL tworząca wskaźnik na obiekt
iteratora do wyszukiwania interpolacji 
\item \verb+gsl_interp_init+ -- metoda GSL wyliczająca funkcję interpolującą
\item \verb+gsl_interp_eval+ -- metoda GSL zwracająca punkt funkcji interpolującej 
\end{itemize}
\medskip{}


Używane zmienne:
\begin{itemize}
\item \verb+double RRmax+ -- najdłuższy interwał RR
\item \verb+double RRmin+ -- najkrótszy interwał RR
\item \verb+IntSignal Histogram+ -- tymczasowy wektor punktów histogramu
\item \verb+unsigned Histogram_size+ -- liczba wszystkich kolumn histogramu
\item \verb+unsigned int RR_intervals_size+ -- liczba wszystkich interwałów RR
\item \verb+double minimum+ -- zmienna przechowująca minimalną wartość całki z algorytmu
wyznaczania TINN
\item \verb+double x[3], y[3]+ -- tablica współrzędnych trójkąta aproksymującego
\end{itemize}


\paragraph{\texttt{GeometricAnalysis::MakePoincareAndSDParams}} -- funkcja wylicza punkty wykresu Poincare oraz parametry: SD1 i SD2.

\medskip{}


Argumenty funkcji:
\begin{itemize}
\item \verb+RR_intervals+ -- wektor interwałów RR w milisekundach
\end{itemize}
\medskip{}


Funkcja zwraca:
\begin{itemize}
\item poincare\_x -- współrzędne osi OX wykresu Poincare
\item poincare\_y -- współrzędne osi OY wykresu Poincare
\item parametr SD1
\item parametr SD2
\end{itemize}
\medskip{}


Używane funkcje:
\begin{itemize}
\item \verb+gsl_vector_get+ -- metoda GSL pobierająca zmiennoprzecinkową wartość
danego elementu z~wektora
\item \verb+gsl_vector_int_set+ -- metoda GSL zapisująca całkowitą wartość w
danym elemencie wektora
\item \verb+gsl_stats_sd+ -- metoda GSL wyliczająca odchylenie standardowe
\item \verb+gsl_vector_int_sub+ -- metoda GSL odejmująca dwa wektory
\end{itemize}
\medskip{}


Używane zmienne:
\begin{itemize}
\item \verb+unsigned int RR_intervals_size+ -- liczba wszystkich interwałów RR
\item \verb+IntSignal diff+ -- wektor różnic pomiędzy sąsiednimi interwałami RR
\end{itemize}
%

\paragraph{\texttt{GeometricAnalysis::SetHRV2Params}} -- funkcja przekazuje wyniki analizy HRV metodami geometrycznymi do instancji
klasy ECGHRV2.

\medskip{}


Argumenty funkcji:
\begin{itemize}
\item \verb+ECGHRV2 &hrv2+
\end{itemize}
\medskip{}


Funkcja zwraca:

Funkcja zwraca wyniki analizy geometrycznej do \verb+ECGHRV2 &hrv2+.

\medskip{}


Używane funkcje:
\begin{itemize}
\item wszystkie funkcje pozwalające na zapis wyników do atrybutów klasy
ECGHRV2. \medskip{}

\end{itemize}
Używane zmienne:

Funkcja nie używa dodatkowych zmiennych.%


\paragraph{\texttt{GeometricAnalysis::MakeRRsignal}} -- funkcja tworząca testowy wektor interwałów RR w milisekundach.

\medskip{}


Argumenty funkcji:

Nie przyjmuje argumentów.

\medskip{}


Funkcja zwraca:

Funkcja zwraca testowy wektor interwałów RR jako instancję klasy ECGRs.

\medskip{}


Używane funkcje:
\begin{itemize}
\item \verb+gsl_vector_int_set+ -- metoda GSL zapisująca całkowitą wartość w
danym elemencie wektora
\item \verb+setRs+ -- zapisuje wektor w instancji klasy ECGRs
\end{itemize}
\medskip{}


Używane zmienne:
\begin{itemize}
\item \verb+int RRsignal_length+ -- długość wektora interwałów RR
\item \verb+int RRsignal[]+ -- tablica wartości interwałów RR
\end{itemize}
%

\paragraph{\texttt{GeometricAnalysis::setParams}} -- funkcja modyfikująca szerokość kolumny histogramu.

\medskip{}


Argumenty funkcji:
\begin{itemize}
\item \verb+ParametersTypes &parameterTypes+ -- zawiera dane ustawień poszczególnych
modułów
\end{itemize}
\medskip{}


Funkcja zwraca:

Nie zwraca niczego. 

\medskip{}


Używane funkcje:
\begin{itemize}
\item \verb+parameterTypes.find(histogram_bin_length)+
- wyszukuje parametr modyfikujący szerokość kolumny histogramu\medskip{}

\end{itemize}
Używane zmienne:
\begin{itemize}
\item \verb+double HistogramParameter+ -- zmienna przechowująca rezultat funkcji \verb+parameterTypes.find+.
\item \verb+double HistogramBinLength+ -- atrybut klasy \verb+GeometricAnalysis+ przechowujący
szerokość kolumny histogramu.
\end{itemize}


\subsubsection{Warunki Testowania}
\label{sec:hrv2:tests}

Podczas procesu implementacji, moduł HRV2 poddawany był dwóm rodzajom
testów:
\begin{itemize}
\item Pierwszy test polegał na porównaniu wyników działania modułu zaimplementowanego
w Matlabie oraz Visual Studio Ultimate 2010. W tym celu jako dane
wejściowe wykorzystano wektor interwałów RR wczytywany z pliku tekstowego
dla Matlaba z tablicy jednowymiarowej dla Visual Studio.
\item Drugi test sprawdzał działanie programu, gdy ten korzystał z testowych
wyników modułu R\_PEAKS.
\end{itemize}

\subsubsection{Wyniki}
\label{sec:hrv2:results}

Rysunki \ref{fig:hrv2_1}, \ref{fig:hrv2_2} oraz \ref{fig:hrv2_3} prezentują wyniki pierwszego testu opisanego w poprzedniej
podsekcji:

\begin{center}
%
\begin{figure}
\begin{centering}
\includegraphics[scale=0.4]{include/hrv2_1}
\par\end{centering}

\caption{Wyniki modułu HRV2 -- Visual Studio Ultimate 2010}
\label{fig:hrv2_1}
\end{figure}

\par\end{center}

\begin{center}
%
\begin{figure}
\begin{centering}
\includegraphics[scale=0.7]{include/hrv2_2}
\par\end{centering}

\caption{Wykres histogramu -- Matlab}
\label{fig:hrv2_2}
\end{figure}

\par\end{center}

\begin{center}
%
\begin{figure}
\begin{centering}
\includegraphics[scale=0.7]{include/hrv2_3}
\par\end{centering}

\caption{Wykres Poincare -- Matlab}
\label{fig:hrv2_3}
\end{figure}

\par\end{center}

\subsection{HRV DFA}
\label{sec:dfa}
Autorzy: Mikołaj Rzepka i Szczepan Czaicki

\subsubsection{Opis zadania}
\label{sec:dfa:desc}

\begin{description}
\item[Temat] Beztrendowa analiza fluktuacji tachogramu (ang. Detrended Fluctuation Analysis)
\item[Opis] Poszukiwanie wciąż nowych rozwiązań oceny ilościowej tachogramu doprowadziło do
opracowania nowej metody –- beztrendowej analizy fluktuacji, która to jest rozszerzeniem
zwyczajnej analizy fluktuacji. Podejście to sprawdziło się w „odkrywaniu” rozległych
i długotrwających korelacji w sygnale EKG. Zasadniczą zaletą tego podejścia jest możliwość
zastosowania algorytmu dla szeregów czasowych, których podstawowe statystyki mają
charakter niestacjonarny (zmieniają się w czasie). Celem projektu jest opracowanie
i implementacja metody wyznaczającej parametry DFA.
\item[Dane] ciąg próbek załamków R z modułu R\_PEAKS
\item[Szukane] moduł programu wyznaczający, wyświetlający parametry analizy DFA w postaci
tabel oraz rysujący wykresy analizy DFA
\end{description}

\subsubsection{Badania literaturowe}
\label{sec:dfa:papers}
Beztrendowa analiza fluktuacyjna polega na określeniu zależności błędu liniowej interpolacji $F(n)$ dla odpowiednio zdefiniowanego szeregu od wielkości okna interpolacji $n$. Metoda ta bierze pod uwagę różnice w lokalnej zawartości danych i może być stosowana dla jednoczesnej analizy całego szeregu danych zawierającego wiele niejednolitych obszarów. W przeciwieństwie do funkcji f(n) analizy fluktuacyjnej, mającej zmienne nachylenie nawet dla n mniejszego od wielkości typowego obszaru stacjonarnego, funkcja beztrendowa F(n) wykazuje nachylenie stałe przez kilka rzędów wielkości. Metoda DFA znakomicie wykazała  różnice między sekwencjami kodującymi a niekodującymi DNA. Użyto więc metody również do wykrywania korelacji w zjawiskach pogodowych, ekonofizyce (np. fluktuacje wartości indeksów giełdowych) oraz do sygnałów fizjologicznych, przede wszystkim zmienności rytmu serca.
Metoda beztrendowej analizy fluktuacji (ang. Detrended Fluctuation Analysis, DFA) polega
na powtórnej analizie sygnału w blokach o długościach kolejno rosnących, aż do określonej
wartości N. W pierwszej kolejności tachogram jest całkowany zgodnie ze wzorem:
\begin{equation}
\sum_{i=1}^{N}({RR(i)-RR_{avg}})
\end{equation}
gdzie $RR(i)$ to $i$-ty interwał, zaś $RR_{avg}$ to średni interwał między kolejnymi uderzeniami (w dzień średnio około 0.7 sekundy).
Następnie tak scałkowany sygnał podzielony zostaje na bloki o określonych długościach n
i w każdym z bloków następuje wyznaczenie trendu lokalnego $y_n(k)$ metodą najmniejszych kwadratów. Skorzystano w tym celu ze wzorów na współczynniki linii trendu $a$ i $b$:
\begin{equation}
 \frac {S \cdot S_{xy} - S_x\cdot S_{y}}{\Delta}
\end{equation}
\begin{equation}
b = \frac {S_{xx} \cdot S_y - S_x\cdot S_{xy}}{\Delta}
\end{equation}
Wykorzystane we wzorach na współczynniki prostej oznaczenia rozpisano poniżej:
\begin{equation}
S = \sum_{i=1}^n 1 = n
\end{equation}
\begin{equation}
S_x = \sum_{i=1}^n x_i
\end{equation}
\begin{equation}
S_y = \sum_{i=1}^n y_i
\end{equation}
\begin{equation}
S_{xx} = \sum_{i=1}^n x_i^2
\end{equation}
\begin{equation}
S_{xy} = \sum_{i=1}^n x_i y_i
\end{equation}
\begin{equation}
S_{yy} = \sum_{i=1}^n y_i^2
\end{equation}
\begin{equation}
\Delta = S \cdot S_{xx} - (S_x)^2
\end{equation}

Proces jest powtarzany dla coraz większych bloków, aż do pełnej długości N . Wówczas
można wyznaczyć zależność fluktuacji od rozmiaru bloku jako
\begin{equation}
F(n) = \sqrt{\frac{1}{N}\sum_{k=1}^{N}(y(k) - y_n(k))^2}
\end{equation}
Zależność ta może zostać umieszczona na wykresie z logarytmiczno-logarytmicznym. Można
zdefiniować współczynnik skalowania
\begin{equation}
\alpha = \lim_{n \to \infty}\frac{\log F(n)}{\log n}
\end{equation}
W zależności od wartości wyliczonego współczynnika możemy wyciągnąć różne wnioski odnoszące się do sygnału wejściowego:
  \begin{itemize}
  \item Dla 0 $<$ $\alpha$ $<$ 0.5 obserwuje się
długozasięgowe antykorelacje (po
krótszym interwale bardziej
prawdopodobne jest wystąpienie
dłuższego), natomiast dla 0.5 $<$ $\alpha$ $<$ 1
dodatnie korelacje długozasięgowe,
  \item Przypadek $\alpha$ = 1 odpowiada szumowi $1/f$
  \item Dla $\alpha$ $\geq$ 1 korelacje istnieją, lecz przestają
mieć charakter potęgowy,
  \item Wartość $\alpha$ = 1.5 odpowiada scałkowanemu
sygnałowi widma białego, czyli błądzeniu
przypadkowemu
  \end{itemize}


\subsubsection{Opis procedur i metod}
\label{sec:dfa:procs}
Główna funkcjonalność modułu znajduje się w klasie \verb|DFAAnalyzer|. Klasa ta,
według przyjętego schematu rozszerza abstrakcyjny moduł analizy beztrendowej analizy fluktuacji tachogramu.
 

Lista i opis najważniejszych funkcji:

\begin{lstlisting}
void DFAAnalyzer::runModule (const ECGRs &ecgRs, ECGHRVDFA &ecghrvdfa)
\end{lstlisting}

Wirtualna funkcja uruchamiająca moduł HRV DFA. Ustawia instancje klas, w których przechowuje dane i uruchamia fukcję \verb|calcDFA| obliczającą podstawowe parametry związane z beztrendową analiza fluktuacji tachogramu.

Argumenty funkcji:
\begin{itemize}
\item \verb+ECGRs &ecgRs+ -- wektor numerów próbek zawierający załamki R
\item \verb+ECGHRVDFA &ecghrvdfa+ -- instancja klasy przechowującej wyniki analizy HRV DFA
\end{itemize}
\medskip{}

\begin{lstlisting}
void DFAAnalyzer::calcDFA()
\end{lstlisting}

Funkcja obliczająca podstawowe parametry związane z beztrendową analiza fluktuacji tachogramu. Najpierw wylicza z sygnału długość interwałów $RR$, następnie całkuje otrzymany sygnał i wywołuje funkcję \verb|calc_all_trends| wyliczającą parametry linii trendów dla kolejnych wielkości okien. Następnie funkcja przelicza wartości fluktuacji w kolejnych oknach, wpisuje je do instancji klasy \verb|ECGHRVDFA| w celu ich późniejszego wyrysowania. Na końcu obliczana jest wartość współczynnika skalowania $\alpha$, również wpisywana do instancji klasy \verb|ECGHRVDFA|.
\medskip{}

\begin{lstlisting}
trend_coefs** DFAAnalyzer::calc_all_trends(double* rr_integrated,
	int length)
\end{lstlisting}

Funkcja wyliczająca parametry linii trendów dla kolejnych wielkości okien sygnału \verb|rr_integrated|. Dzieli sygnał na okna o wielkości od 2 próbek do całej długości sygnału i dla każdego okna wywołuje funkcję \verb|calc_trend| obliczającą parametry $a$ i $b$ linii trendu metodą najmniejszych kwadratów.

Argumenty funkcji:
\begin{itemize}
\item \verb+double* rr_integrated+ -- sygnał (tablica) interwałów $RR$
\item \verb+int length+ -- długość tablicy interwałów $RR$
\end{itemize}
\medskip{}

\begin{lstlisting}
trend_coefs* DFAAnalyzer::calc_trend(double* rr_integrated, int start,
	int w_length)
\end{lstlisting}

Funkcja obliczająca parametry $a$ i $b$ linii trendu metodą najmniejszych kwadratów dla podanego sygnału od próbki \verb|start| w oknie długości \verb|w_length|.

Argumenty funkcji:
\begin{itemize}
\item \verb+double* rr_integrated+ -- sygnał (tablica) interwałów $RR$
\item \verb+int start+ -- numer próbki sygnału interwałów $RR$ od którego funkcja zaczyna przeliczenia
\item \verb+int w_length+ -- ilość próbek sygnału interwałów $RR$ brana do przeliczeń począwszy od próbki \verb|start|
\end{itemize}
\medskip{}

\subsubsection{Warunki Testowania}
\label{sec:dfa:tests}

Grupa wykonująca zadania były uzależniona od wyników grupy implementującej moduł R\_PEAKS. W związku z tym, aby uniezależnić się od pracy tej grupy, do testowania algorytmu używano sygnałów zawierających piki R pobranych z internetu. Wyniki porównywano z wynikami uzyskiwanymi z~wyliczeń skryptu napisanego w Matlabie w poprzednim semestrze.

\subsubsection{Wyniki}
\label{sec:dfa:results}

W pierwszym etapie prac porównywano wyniki pracy algorytmu do wyników uzyskiwanych z wyliczeń skryptu napisanego w Matlabie, udało się osiągnąć takie same rezultaty. Na efekty wizualizacji wyników czekaliśmy do samego końca, z powodu braku czasu zaimplementowano wizualizację jednego wykresu: zależności logarytmu fluktuacji od logarytmu z wielkości okna.

\subsection{Klasyfikacja zespołów QRS}
\label{sec:qrs_class}
Autorzy: Krzysztof Bębenek

\subsubsection{Opis zadania}
\label{sec:qrs_class:desc}

\begin{description}
\item[Temat] Metody detekcji morfologicznego pochodzenia zespołu QRS.
\item[Opis] Proces automatycznego klasyfikowania zespołów QRS należy do jednego z trudniejszych procesów podczas przetwarzania sygnałów EKG. Jest jednak on niezbędny w dalszych etapach analizy podczas których brane są pod uwagę tylko niepoprawne pobudzenia. Do prostych sposób oceny morfologi pobudzenia służą:
  \begin{itemize}
  \item współczynnik kształtu Malinowskiej,
  \item stosunek części ujemnej do dodatniej sygnału,
  \end{itemize}
\item[Dane] ciąg próbek przefiltrowanego sygnału EKG, wektory QRS\_onset oraz QRS\_end określające początek oraz koniec zespołu QRS
\item[Szukane] moduł programu klasyfikujący zespołu QRS na podstawie ich morfologii 
\end{description}

\subsubsection{Badania literaturowe}
\label{sec:qrs_class:papers}

\subsubsection{Opis procedur i metod}
\label{sec:qrs_class:procs}

\subsubsection{Warunki testowania}
\label{sec:qrs_class:tests}

\subsubsection{Wyniki}
\label{sec:qrs_class:results}



\subsection{ST interval}
\label{sec:st_interval}
Autorzy: Bartłomiej Bułat i Krzysztof Piekutowski

\subsubsection{Opis zadania}
\label{sec:st_interval:desc}
Celem projektu jest wyznaczenie odcinków ST, pomiar poziomu odcinka ST względem
linii izoelektrycznej i jego nachylenia, a także detekcja epizodów ST oraz
parametry ilościowe i jakościowe epizodów ST. Celowość analizy odcinka ST
względem linii izoelektrycznej związana jest z predykcją choroby wieńcowej,
miażdżycy oraz niedotlenieniem mięśnia serca.

Dane przyjmowane przez moduł:
\begin{itemize}
  \item sygnał ECG\_BASELINE;
  \item wektor numerów próbek załamków R z modułu R PEAKS;
  \item wektor numerów próbek punktów charakterystycznych z modułu WAVES:
  \begin{itemize}
    \item $QRS_{onset}$
    \item $QRS_{end}$
    \item $T_{end}$
  \end{itemize}  
\end{itemize}

Dane zwracane przez moduł:
\begin{itemize}
  \item wektor wykrytych odcinków ST zawierających informacje:
  \begin{itemize}
    \item początek i koniec odcinka;
    \item poziom izolinii;
    \item pomiar poziomu względem izolinii;
    \item pomiar nachylenia odcinka ST;
    \item opis słowny interwału;
  \end{itemize}
  \item wektor wykrytych epizodów wraz z parametrami:
  \begin{itemize}
    \item numer próbki początku epizodu;
    \item numer próbki końca epizodu.
  \end{itemize}
\end{itemize}

\subsubsection{Badania literaturowe}
\label{sec:st_interval:papers}

Analiza odcinka ST w sygnale EKG jest kluczowa we wczesnym wykrywaniu chorób
wieńcowych, chorobach niedokrwiennych i miażdżycy. Badając uniesienie i
nachylenie odcinka względem izolini można uzyskać informacje na temat
zagrożenia wymienionymi chorobami i podjąć wczesną interwencję. Dość sporym
mankamentem tych parametrów diagnostycznych jest ich charakter amplitudowy,
gdyż wiadomo, że wielkość amplitudy zależy często od wielu pozakardiologicznych
przyczyn, m.in. jakości kontaktu elektrod ze skórą.

Książka \cite{AUGUST1} prezentuje bardzo prosty algorytm diagnostyki odcinka
ST. Algorytm składa się z~trzech kroków: lokalizacji odcina, obliczenia uniesienia i
nachylenia oraz wykrywanie epizodów.

\begin{figure}[h!]
  \centering
  \includegraphics[width=0.6\textwidth]{include/ecg_st}
  \caption{Kluczowe punkty analizy odcinka ST. ISO -- punk referencyjny leżący
  na izolinii, J -- początek odcinka, ST -- koniec odcinka. Opracowanie własne.}
  \label{fig:ecg_st}
\end{figure}

Na rysunku \ref{fig:ecg_st} przedstawiono kluczowe punkty, które należy
wyznaczyć przed wyznaczaniem parametrów diagnostycznych odcinka ST. Według
wyżej wymienionego algorytmu, punktem referencyjnym w założeniu znajdującym się
na izolinii jest początek zespołu QRS i koniec odcinka P. Punkt J, który jest
początkiem badanego odcinka i znajduje się $45ms$ za szczytem załamka R. Punkt
ST będący końcem badanego odcinka znajduje się $60ms$ później. Tych stałych
wartości używa się w prostych algorytmach, dodatkowo można uwzględnić zmienność
rytmu serca i uzależnić przesunięcie tych wartości od aktualnego odstępu RR.
Uniesieniem lub obniżeniem odcinka jest różnica wartości miedzy punktem ISO, a
punktem ST. Nachyleniem prostej przechodzącej przez punkty J i ST nazywamy
nachyleniem odcinka. Prawidłowy odcinek ST powinien być poziomy i nie być
uniesiony ani obniżony. Algorytm wykrywania epizodów polega na wykryciu takich
okresów zapisu EKG w którym przez czas dłuższy niż 60 sekund występują ciągle
nieprawidłowe odcinki ST.

W artykule \cite{SHEN1} przedstawiono bardziej złożoną procedurę wyznaczania
punktów charakterystycznych odcinka. W oparciu o znalezione wcześniej punkty
ISO, J i punkt końca odcinka T ($T_{end}$) znajduje szczyt załamka T z
wykorzystaniem dekompozycji falkowej. następnie wyznacza punkt ST jako
najbardziej odległą wartość sygnału od prostej łączącej punkt J i punkt
szczytowy załamka T.

Obniżenie załamka jest wyznaczane metodą J+X, który polega na wyznaczeniu,
podobnie jak w~pierwszym algorytmie, różnicy miedzy punktem J+X (punkt X
milisekund za punktem J), a punktem ISO. X jest wyznaczany na podstawie
aktualnego HRV. Ocena przesunięcia odcina ST względem izolinii odbywa się z
użyciem progu K, dla bazy MIT-BIH użyta wartością może być 0.01. Odcinek uznany
jest za podniesiony gdy przesunięcie jest większe od K, obniżony gdy jest
mniejszy od -K i~normalny, prawidłowy w pozostałych przypadkach.

Następnie w zalezności od oceny przesunięcia odcinka wyznaczany jest
referencyjny punkt końcowy. Dla odcinka uniesionego jest to punkt szczytowy
fali T, dla normalnego i obniżonego -- punkt ST. Kolejno, na podstawie
charakteru odcinka do J20 (punkt J przesunięty o 20ms) do końcowego punktu
referencyjnego okreslane jest czy odcinek jest prosty czy zakrzywiony, a
następnie określany jest charakter zakrzywienia (wklęsły/wypukły) lub charakter
monotoniczny odcinka prostego.

Rodzaj zakrzywienia odcinka ST wyznaczany jest na podstawie stosunku liczby
punktów badanego fragmentu sygnału które są ponad/poniżej liniowej interpolacji
odcinka do liczby wszystkich punktów odcinka. Jeżeli ten stosunek dla punków
ponad prostą jest większy od 70\% to odcinek jest wypukły, jeżeli 70\% punktów
jest poniżej prostej, odcinek jest wklęsły.

Charakter monotoniczny prostego odcinka ST jest określany na podstawie
nachylenia interpolacji liniowej badanego fragmentu sygnału.

Epizody ST są wyznaczane podobnie jak w pierwszym algorytmie.

Podczas implementacji wprowadzono kilka poprawek, z których należy
wyszczególnić:

Zmiana sposobu wyznaczania szczytu fali T. Zamiast dekompozycji falkowej użyto
pochodnej sygnału. Ostatni punkt przejścia pochodnej przez zero jest szczytem
fali.

Porzucono branie pod uwagę HRV, przy wyznaczaniu kolejnych punktów.
\subsubsection{Opis procedur i metod}
\label{sec:st_interval:procs}

Główna funkcjonalność modułu znajduje się w klasie \verb|STAnalysis|. Klasa ta,
według przyjętego schematu rozszerza abstrakcyjny moduł analizy odcinka ST. W
tej klasie znajdują się drzewo klas prywatnych reprezentujących poszczególne
algorytmy analizy. W szczególności dwie klasy \verb|SimpleAnalizer| oraz
\verb|ComplexAnalizer|. Pierwsza klasa reprezentuje najprostszy algorytm
zaprezentowany w \cite[p.~155]{AUGUST1}, druga zaś implementuje zmodyfikowany
algorytm opisany w \cite{SHEN1}.

%% TODO: mienić nazwę funkcji w programie
Lista i opis najważniejszych funkcji:

\begin{lstlisting}
void STAnalysis::SimpleAnalizer::analyse(const int it, const ECGRs& rpeaks,
  const ECGWaves& waves, const ECGSignalChannel& signal,
  const  ECGInfo& info, ECGST& output);
\end{lstlisting}

Funkcja analizująca załamek ST w zespole QRS numer \verb|it|. Wykorzystując
punkty sygnału obliczone we wcześniejszych modułach (tj. $QRS_{onset}$,
$QRS_{end}$ oraz $T_{end}$) oblicza koniec odcinka ST i wylicza jego parametry:
\begin{itemize}
  \item przesunięcie względem izolinii
  \item wartość nachylenia w stosunku do izolini
  \item długość odcinka ST
\end{itemize}
Używając wcześniej ustawionego parametru \verb|simple_thresh|, odcinek ST
określany jest jako ,,normalny'', ,,uniesiony'' lub ,,obniżony''.

Parametry:
\begin{itemize}
  \item \verb|const int it| -- numer badanego zespołu QRS, jako liczba
    porządkowa numerów próbek z tablicy załamków R.
  \item \verb|const ECGR& rpeaks| -- struktura zawierająca tablicę numerów
    próbek kolejnych załamków R
  \item \verb|const ECGEaves& waves| -- struktura zawierająca tablice
    przechowujące numery próbek punktów charakterystycznych kolejnych zespołów
    QRS: $P_{onset}$, $QRS_{onset}$, $QRS_{end}$ oraz $T_{end}$.
  \item \verb|const ECGSignalChanel& signal| -- jeden kanał odfiltrowanego
    sygnału z usuniętym przesunięciem izolinii.
  \item \verb|const ECGInfo& info| -- struktura zawierająca informacje o badanym
    sygnale EKG, m.in. częstotliwość.
  \item \verb|ECGST& output| -- parametr wyjściowy, struktura zawierająca
    tablice wszystkich interwałów ST wraz z ich parametrami, oraz tablice
    zarejestrowanych epizodów ST wraz z ich parametrami.
\end{itemize}

\begin{lstlisting}
void STAnalysis::ComplexAnalizer::analyse(const int it, const ECGRs& rpeaks,
  const ECGWaves& waves, const ECGSignalChannel& signal,
  const  ECGInfo& info, ECGST& output);
\end{lstlisting}

Funkcja analizująca załamek ST w zespole QRS numer \verb|it|. Wykorzystując
punkty sygnału obliczone we wcześniejszych modułach (tj. $QRS_{onset}$,
$QRS_{end}$ oraz $T_{end}$) oblicza koniec odcinka ST (wykorzystując bardziej
zaawansowane algorytmy od poprzedniej funkcji), a następnie oblicza jego parametry:
\begin{itemize}
  \item przesunięcie względem izolinii,
  \item wartość nachylenia w stosunku do izolinii,
  \item długość odcinka ST,
  \item klasyfikacja kształtu.
\end{itemize}
Używając wcześniej ustawionego parametru \verb|complex_thresh|, odcinek ST
określany jest jako ,,normalny'', ,,uniesiony'' lub ,,obniżony''. Algorytm
ocenia również to, czy załamek jest prosty czy zakrzywiony w oparciu o
wcześniej ustawiony parametr \verb|type_thresh|. Dla prostych odcinków określa,
czy kierunek narasta, opada czy jest poziomy. Do określenia tej cechy
wykorzystywany jest parametr \verb|slope_thresh|. Dla zakrzywionych odcinków ST
oceniana jest wypukłość krzywej. Szczegółowy opis działania znajduje się w
poprzednim rozdziale.

Parametry:
\begin{itemize}
  \item \verb|const int it| -- numer badanego zespołu QRS, jako liczba
    porządkowa numerów próbek z tablicy załamków R.
  \item \verb|const ECGR& rpeaks| -- struktura zawierająca tablicę numerów
    próbek kolejnych załamków R
  \item \verb|const ECGEaves& waves| -- struktura zawierające tablice
    przechowujące numery próbek punktów charakterystycznych kolejnych zespołów
    QRS: $P_{onset}$, $QRS_{onset}$, $QRS_{end}$ oraz $T_{end}$.
  \item \verb|const ECGSignalChanel& signal| -- jeden kanał odfiltrowanego
    sygnału z usuniętym przesunięciem izolinii.
  \item \verb|const ECGInfo& info| -- struktura zawierająca informacje o badanym
    sygnale EKG, m.in. częstotliwość.
  \item \verb|ECGST& output| -- parametr wyjściowy, struktura zawierająca
    tablice wszystkich interwałów ST wraz z ich parametrami, oraz tablice
    zarejestrowanych epizodów ST wraz z ich parametrami.
\end{itemize}

Funkcje pomocnicze klasy \verb|ComplexAnalizer|:

\begin{lstlisting}
int STAnalysis::ComplexAnalizer::getTPeak(const OtherSignal& sig,
  int from, int to);
\end{lstlisting}

Funkcja wyszukująca położenia punktu $T_{peak}$ na zadanym odcinku sygnału.
Początkiem wyszukiwania szczytu fali T jest zwykle punkt 20ms za punktem
$QRS_{end}$, a punktem końcowym jest koniec fali T. Funkcja wykorzystuje
dyskretną pochodną sygnału.

Parametry:
\begin{itemize}
  \item \verb|const OtherSignal& sig|  -- cały, odfiltrowany sygnał EKG
  \item \verb|int from|  -- numer próbki, w której należy zacząć poszukiwania
  \item \verb|int to|  -- numer próbki, w której należy skończyć poszukiwania
\end{itemize}


\begin{lstlisting}
std::pair<int, double> maxDistanceSample(const OtherSignal& signal, 
  int from, int to);
\end{lstlisting}

Funkcja szukająca numeru próbki najbardziej oddalonego od liniowej interpolacji
fragmentu sygnału miedzy dwoma punktami. Oprócz numeru próbki od początku
sygnału, zwracana jest wartość największej różnicy miedzy punktem sygnału, a
prostą interpolacji.

Parametry:
\begin{itemize}
  \item \verb|const OtherSignal& signal|  -- badany sygnał
  \item \verb|int from|  -- numer próbki będącej początkiem interesującego
    fragmentu sygnału, równocześnie, pierwszy punkt interpolacji liniowej.
  \item \verb|int to|  -- numer próbki będącej końcem interesującego fragmentu
    sygnału, równocześnie, drugi punkt interpolacji
\end{itemize}

\begin{lstlisting}
std::pair<int, int> overBelowSamples(const OtherSignal& signal,
  int from, int to);
\end{lstlisting}

Funkcja obliczająca ilość próbek ponad i poniżej prostej interpolującej sygnał
miedzy dwoma punktami. Funkcja wykorzystywana jest do określenia wypukłości
odcinka ST.

Parametry:
\begin{itemize}
  \item \verb|const OtherSignal& signal|  -- badany sygnał
  \item \verb|int from|  -- numer próbki będącej początkiem interesującego
    fragmentu sygnału, równocześnie, pierwszy punkt interpolacji liniowej.
  \item \verb|int to|  -- numer próbki będącej końcem interesującego fragmentu
    sygnału, równocześnie, drugi punkt interpolacji
\end{itemize}

\subsubsection{Warunki testowania}
\label{sec:st_interval:tests}
W trakcie implementacji do testowania algorytmów używaliśmy surowego sygnału MIT BIH oraz wzorcowych załamków R z pliku z rozszerzeniem atr dostępnego razem z sygnałem. Testy polegały na wizualnej ocenie poprawności wykrytych punków przedstawionych na wykresie. Kolejnym krokiem było testowanie algorytmu na przefiltrowanym sygnale oraz załamkach R zwracanych przez moduł RPeaks, również na bazie oceny wizualnej wyznaczonych właściwości odcinka ST. Flaga \verb|DEVELOPMENT| pozwala na przełączanie się między użyciem załamków wzorcowych (gdy jest ustawiona) i wykrywanych przez moduł RPeaks (gdy nie jest ustawiona).

\subsubsection{Wyniki}
\label{sec:st_interval:results}

\begin{figure}[h!]
  \centering
  \includegraphics[width=0.9\textwidth]{include/st_gui}
  \caption{Zakładka modułu w GUI przedstawiająca przykładowy zespół QRS z
  zaznaczonymi punktami charakterystycznymi analizy odcinka ST. }
  \label{fig:st_gui}
\end{figure}

%%TODO można tu dodać printscreen wykresu ST, gdy zacznie działać aplikacja
Resultatem projektu są dwa algorytmy wyznaczania i analizy parametrów odcinków ST. Prostrzy algorytm opisany w \cite{AUGUST1} został w całości zaimplementowany i przetestowany. Drugi, rozszerzony algorytm z \cite{SHEN1} również został zrealizowany, jednak ze względu na złożoność matematyczną niezbędne okazało się uproszczenie jednego z kroków. Algorytm ten zakłada wykorzystanie 5-warstwowej dekompozycji falkowej splinów 4-stopnia do detekcji punktu $T_{peak}$, co okazało się być niewykonalne technicznie. Wykrycie tego punktu zostało zrealizowane przy użyciu filtru różnicowego wyszukującego maximum lokalne sygnału. Poprawność wyników z obu algorytmów w dużym stopniu zależy od jakości wyników działania modułów nadrzędnych, dostarczających przefiltrowany sygnał EKG oraz zespoły QRS.
Na rysunku \ref{fig:st_gui} przedstawiono przykładowy wynik modułu analizy ST.


\subsection{T wave alt}
\label{sec:t_wave_alt}

\subsubsection{Opis zadania}
\label{sec::t_wave_alt:desc}

\begin{description}
\item[Temat] Ocena ilościowa alternansu załamka T
\item[Opis] W ostatnich latach alternans załamka T stał się
 elementem szczególnych zainteresowań klinistów w środowisku kardiologów.
 Dowiedziono, że pojawienie się naprzemienności kształtu, amplitudy czy
 polarności załamka T jest zapowiedzią zaburzeń komorowych rytmu serca 
(arytmii komorowych). Wykazano również korelację amplitudy alternansu 
załamka T i spadku progu trzepotania przedsionków i komór. Celem projektu
 jest oznaczenie załamków T, w których wykryto alternans oraz ich ilościowa analiza.
\item[Dane] ciąg punktów charakterystycznych z modułu WAVES i sygnał 
\item[Szukane] moduł programu oznaczający załamki T, w których wykryto alternans
 oraz dokonujący ich ilościowej analizy. Wszystkie dane należy odpowiednio nanieść
 na wykresy i przedstawić w~postaci graficznej (parametry alternansu).
\end{description}

\subsubsection{Badania literaturowe}
\label{sec::t_wave_alt:papers}
Alternans załamka T to naprzemienna zmiana wektora oraz amplitudy załamka T w sygnale EKG. Zmiany wielkości są mierzone w mikrowoltach, dlatego do analizy sygnału i detekcji TWA wykorzystywany jest sprzęt o bardzo dużej czułości. 

Detekcja alternansu załamka T jest obecnie przedmiotem wielu badań oraz analiz prowadzonych przez ośrodki kardiologiczne, laboratoria w ośrodkach akademickich oraz centra badawcze na całym świecie. Nie istnieje obecnie jedna uniwersalna metoda dająca obiektywne wyniki i mogącą z 100\% pewnością stwierdzić, że dany przebieg EKG zawiera alternans. Dodatkowo testy skuteczności różnych metod bardzo często dają sprzeczne ze sobą wyniki. W związku z faktem istnienia bardzo dużej różnorodności metod, narzędzi oraz algorytmów mających na celu wykrycie alternansu załamka T postanowiono przeprowadzić badania umożliwiające ocenę różnych metod.

Na podstawie przeanalizowanej literatury oraz algorytmów biorących udział w konkursie \verb|Physionet|  \verb|Challenge 2008|, zostały wyodrębnione trzy główne metody, które zostały ostatecznie wzięte pod uwagę przy wyborze algorytmu detekcji oraz projektowaniu modułu:
\begin{itemize}
	\item Moving four-beat window techniqure
	\item Fast Fourier Transform
	\item Spectral analysis PCA
\end{itemize}

Po przeprowadzeniu testów oraz analizie wniosków z pracy  \verb|T-Wave Alternans:|  \verb|A Comparison of|  \verb|Different Measurement Techniques| zdecydowano wybrać podejście  \verb|Moving four-beat window| jako najbardziej odpowiednie do stawianego przed autorami problemu i dającego względnie najlepsze wyniki. O złożoności problemu detekcji TWA świadczy między innymi fakt, że podczas przeprowadzania testów trzech wyżej wymienionych algorytmów jedynie 7 razy wszystkie metody były zgodne co do występowania w sygnale alternansu załamka T. 

Algorytm \verb|Moving four-beat window techniqure| analizuje jednocześnie 4 występujące obok siebie załamki T (okno o rozmiarze 4). W przypadku, gdy zajdą w takim oknie następujące zależności: 
\begin{itemize}
	\item A2-A4 \textless A1-A2
	\item A2-A4 \textless A3-A4
	\item A2-A4 \textless A1-A2
	\item A1-A3 \textless A3-A4
\end{itemize}

dane okno brane jest pod uwagę jako posiadające alternans. Okno jest przesuwane, aż do napotkania końca sygnału. W przypadku gdy 5\% wszystkich okien spełnia wyżej wymieniony warunek, dany sygnał kwalifikowany jest jako zawierający alternans załamka T. 

\begin{center}
%
\begin{figure}
\begin{centering}
\includegraphics[scale=0.4]{include/twa_1}
\par\end{centering}

\caption{Algorytm Movign four-beat window techniqure}
\label{fig:twa_1}
\end{figure}

\par\end{center}

\subsubsection{Opis procedur i metod}
\label{sec::t_wave_alt:procs}

Pełna funkcjonalność modułu została za implementowana w klasie \verb|TWaveAltDetector|. Klasa ta dziedziczy po abstrakcyjnej klasie  \verb|TWaveAltModule|. Moduł uruchamiany jest za pomocą funkcji: 

\paragraph{\texttt{TWaveAltDetector::runModule(const ECGWaves, const ECGSignalChannel, ECGInfo, ECGTWave)}}

\medskip{}

Argumenty funkcji:

\begin{itemize}
\item \verb+ECGWaves ecgWaves+ -- instancja klasy przechowującej wyniki detekcji punktów charakterystycznych sygnału EKG takich jak QRS-onset, QRS-end, T-end, P-onset, P-end;
\item \verb+ECGSignalChannel ecgSignal+ -- instancja klasy zawierajacej badany sygnał EKG;
\item \verb+ECGInfo ecgInfo+ -- informacje na temat wczytanego sygnału;
\item \verb+ECGTWave ecgTWave+ -- iinstancja klasy przechowującej wyniki detekcji oraz wektor numerów próbek sygnału zawierających alternans załamka T;
\end{itemize}
\medskip{}

Używane funkcje:
\begin{itemize}
\item \verb+detectTWaveAlt()+ --główna funkcja klasy odpowiedzialna za całą logikę modułu;
\end{itemize}

\medskip{}

Używane zmienne:
\begin{itemize}
\item \verb+filteredSignal+ --  badany sygnał
\item \verb+wavesDatal+ -- wektor numerów próbek zawierających QRS-onset, QRS-end, T-end, P-onset, P-end
\item \verb+tWaveAltData+ -- wektor numerów próbek zawierających alternans
\end{itemize}


\medskip{}

Detekcja oraz obliczenie parametrów oceny ilościowej alternansu wykonywane są w metodzie  \verb|detectTWaveAlt()|. Cały algorytm przetwarzania sygnału, detekcji alternansu (zgodnie z algorytmem  \verb|Moving four-beat window techniqure|) oraz obliczenia ilościowe związane z analizowanym sygnałem zostały zawarte właśnie w tej funkcji. Funkcja nie przyjmuje żadnych parametrów wejściowych. Funkcja zwraca wartość logiczną 1 w przypadku wykrycia w sygnale alternansu oraz wartość 0 w~przypadku braku alternansu w sygnale. Funkcja oblicza następujące parametry ilościowe alternansu:
\begin{itemize}
\item \verb+ilość okien w sygnale zawierających alternans+
\item \verb+procentowy udział wszystkich okien o rozmiarze 4 zawierających alternans w sygnale+
\end{itemize}


\paragraph{\texttt{bool TWaveAltDetector::detectTWaveAlt()}}

\medskip{}

Argumenty funkcji:

\begin{itemize}
\item \verb+brak+
\end{itemize}
\medskip{}

Używane funkcje:
\begin{itemize}
\item \verb+gsl_vector_int_get()+ --funkcja GSL pobierająca całkowitą wartość danego elementu z wektora
\item \verb+gsl_vector_int_alloc()+ --funkcja GSL tworząca wskaźnik do nowo utworzonego obiektu 
\item \verb+gsl_vector_int_set()+ --funkcja GSL zapisująca całkowitą wartość danego elementu wektora
\item \verb+GetT_end()+ --funkcja zwracająca wektor numerów próbek zawierających załamek T
\item \verb+setTWaveAlt()+ --funkcja ustawiająca wektor numerów próbek zawierających alternans załamka T
\end{itemize}


\medskip{}

Używane zmienne:
\begin{itemize}
\item \verb+filteredSignal+ --  badany sygnał
\item \verb+wavesDatal+ -- wektor numerów próbek zawierających QRS-onset, QRS-end, T-end, P-onset, P-end
\item \verb+tWaveAltData+ -- wektor numerów próbek zawierających alternans
\end{itemize}


\medskip{}

\subsubsection{Warunki testowania}
\label{sec::t_wave_alt:tests}
Wszystkie algorytmy brane pod uwagę przy wyborze docelowego rozwiązania zostały wstępnie zaimplementowane w środowisku Matlab. Po wybraniu odpowiedniego i najbardziej skutecznego rozwiązania (\verb|Moving four-beat window techniqure|) moduł \verb|TWaveAlt| wraz z detekcją załamków T został zaimplementowany w Matlab, a następnie w C++. 

Zostały przeprowadzone następujące testy:

\begin{itemize}
\item Testy różnych algorytmów detekcji załamka T: Zostały opracowane prototypy różnych alogtytmów detekcji alternansu w Matlab. 
\item Testy poprawności działania wybranego algorytmu: Porównanie działania modułu zimplementowanego w Matlab z algorytmem w Visual Studio. 
\item Testy skuteczności wykrywania alternansu T: Zostały przetestowanie sygnały zawierające alternans T, które przedmiotem badań podczas konkursu "Physionet Challenge 2008". Uzyskane wyniki zostały porównane z oficjalnymi rezultatami prac konkursowych.
\end{itemize}

\subsubsection{Wyniki}
\label{sec::t_wave_alt:results}
Przeprowadzone rysunki przedstawiają testy różnych algorytmów detekcji załamka T.

\begin{center}
%
\begin{figure}
\begin{centering}
\includegraphics[scale=0.4]{include/twa_2}
\par\end{centering}

\caption{Test 1  -- TWA algorytmy}
\label{fig:twa_2}
\end{figure}

\par\end{center}


\begin{center}
%
\begin{figure}
\begin{centering}
\includegraphics[scale=0.4]{include/twa_3}
\par\end{centering}

\caption{Test 2  -- TWA skuteczność}
\label{fig:twa_3}
\end{figure}

\par\end{center}

\subsection{HRT}
\label{sec:hrt}

\subsubsection{Opis zadania}
\label{sec:hrt:desc}

Turbulencja rytmu serca jest nową, jeszcze nie do końca poznaną techniką służącą do oceny odpowiedzi węzła zatokowego na przedwczesny skurcz komorowy. Podejście okazało się trafnym narzędziem predykcji nagłych zgonów pacjenta.
Celem projektu była implementacja modułu programu wyznaczającego podstawowe parametry turbulencji rytmu serca (m. in. początek turbulencji, nachylenie turbulencji, itp…).
\\
\\
Autorzy: Łukasz Kutrzuba, Mateusz Krasucki
\\
\subsubsection{Badania literaturowe}
\label{sec:hrt:papers}

\textbf{Przedwczesne pobudzenie komorowe}\\
Przedwczesny skurcz komorowy (ang. Premature Ventricular
Contraction) powstaje na skutek przedwczesnego rozprzestrzenienia się
impulsu w komorach powodującego skurcz mięśnia komór (impuls w~normalnym przypadku powinien pochodzić z węzła zatokowo-przedsionkowego).
Jest to jedna z najczęstszych arytmii.
\textbf{Turbulencja rytmu serca}\\
Turbulencja rytmu serca (ang. Heart Rate Turbulence) jest to reakcja fizjologiczna będąca powrotem do stanu równowagi po przedwczesnym skurczu
komorowym. Analiza sygnałów HRT pozwala na ocenienie jak
szybko i jak energicznie serce reaguje na przedwczesny skurcz
komorowy. Duże znaczenie ma fakt, że występowanie tego
ostatniego zjawiska fizjologicznego jest naturalne dla ogółu osób
dorosłych, co pozwala na wykorzystanie oceny HRT dla
diagnostyki stanu zdrowia. Analiza HRT jest zwykle wykonywana
na długich, 24 godzinnych zapisach EKG.
\par W zapisie reakcji organizmu na PVC można wskazać charakterystyczne etapy. Po przedwczesnym skurczu komorowym następuje przerwa wyrównawcza, następnie faza akceleracji, deceleracja i ostatecznie powrót do normalnego rytmu.
Nie zdołano dotychczas określić przyczyn zjawiska turbulencji rytmu serca, istnieją przypuszczenia, że powoduje je odruch z baroreceptorów.
\par Wystąpienie turbulencji rytmu serca świadczy o prawidłowej reakcji 
układu przywspółczulnego i~zdolności serca do reagowania na przedwczesne pobudzenie -- z tego powodu wystąpienie skróconych interwałów po VPC i ich stopniowe wydłużanie jest uznawane za prawidłową reakcję. Z kolei pojawienie się interwałów dłuższych niż przed 
przedwczesnym pobudzeniem jest negatywnym prognostykiem, świadczy o braku
reakcji organów kontrolujących pracę serca.

\begin{figure}[h!]
\centering
%\includegraphics[width=13cm]{include/HRTSample.eps}
\caption{Wykres przedstawiający wykresy HRT dla zdrowego i chorego człowieka. 
\textit{Źródło: http://www.h-r-t.org, Working Group of Biological Signal Analyses, Technische Universität München Źródło: Axel Bauer, Marek Malik, Georg Schmidt: „Heart Rate Turbulence: Standards of Measurement, Physiological Interpretation, and Clinical Use”, Journal of the American College of Cardiology 2008, tom 52, nr 17} }
\end{figure}

\textbf{Znaczenie diagnostyczne turbulencji rytmu serca}\\
Analiza turbulencji rytmu serca jest jedną z najważniejszych metod pozwalającą na ocenę ryzyka zgonu pacjentów po przebytym zawale, obok frakcji wyrzutowej lewej komory, ilości przedwczesnych skurczów komorowych oraz analizy rytmu serca. Dowiedziono wysokiej wartości predykcyjnej tej metody badania.
Przy analizie wyników pacjentów zostaje dokonane przyporządkowanie do jednej trzech grup -- oba parametry poprawne, jeden parametr negatywny, oba parametry negatywne.
Na przykład w przypadku obu parametrów negatywnych ryzyko zgonu w ciągu roku od przebytego zawału wynosi powyżej 32 procent.

\textbf{Parametry turbulencji rytmu serca}\\
W przypadku turbulencji rytmu serca wykorzystuje się zwykle dwa parametry
-- początek turbulencji (turbulence onset) oraz nachylenie turbulencji (turbulence slope). Turbulence Onset (TO) określa procentową różnicę między częstością serca przedwczesnym skurczu komorowym i częstością serca przed tym skurczem. Obliczanie TO opiera się na różnicy pomiędzy sumą długości dwóch interwałów RR występujących bezpośrednio po pauzie wyrównawczej a sumą dwóch interwałów występujących przed pauzą wyrównawczą.
Wartości dodatnie TO oznaczają zwolnienie, natomiast ujemne przyspieszenie rytmu zatokowego. Pozytywnym prognostykiem są wartości ujemne.

\begin{figure}[h!]
\centering
%\includegraphics[width=15cm]{include/to1.eps}
\caption{Przykładowy schemat ilustrujący wyznaczanie parametru Turbulence Onset.
\textit{ Źródło: http://www.h-r-t.org, Working Group of Biological Signal Analyses, Technische Universität München} }
\end{figure}
\par Turbulence Slope (TS) odpowiada najbardziej stromej linii regresji przeprowadzona przez każde 5 kolejnych odstępów prawidłowego rytmu w tachogramie (spośród 20 kolejnych interwałów występujących zaraz po przerwie wyrónawczej). Obliczenia Turbulence Slope opierają się na tachogramie uśrednionym z wszystkich fragmentów rytmu zatokowego po skurczu komorowym. Warości TS wyrażone w ms/odstęp RR. Prawidłowa wartość parametru to 2.5ms/RR.
\begin{figure}[h!]
\centering
%\includegraphics[width=15cm]{include/ts.eps}
\caption{Przykładowy schemat ilustrujący wyznaczanie parametru Turbulence Slope.
\textit{ Źródło: http://www.h-r-t.org, Working Group of Biological Signal Analyses, Technische Universität München} }
\end{figure}
\newpage

\subsubsection{Koncepcja rozwiązania}
\label{sec:hrt:procs}

\textbf{Detekcja PVC}\\
Przy wyliczeniu parametrów TO i TS brane są pod uwagę niektóre z fragmentów zapisu interwałów RR  -- zwykle jest to fragment zapisu liczący 27 kolejnych 
interwałów ( 5 interwałów przed VPC, interwał z~VPC, pauza kompensacyjna,
około 20 interwałów po pauzie). Dodatkowo, wybierane fragmenty zapisu muszą  charakteryzować się się brakiem występowania nierównego rytmu serca przed i po PVC oraz brakiem zbyt krótkich oraz zbyt długich interwałów RR.
\par W celu selekcji odpowiednich odcinków zastosowano następujący algorytm:
\begin{itemize}
\item Pobranie sygnału z długościami poszczególnych interwałów
\item Każdy kolejny interwał RR jest traktowany jako
potencjalny interwał VPC (z kolei 5 interwałów poprzedzających oraz 20 interwałów występujących po tym interwale będzie wykorzystanych do oceny
sygnału jako poprawny VPC). 
W przypadku, jeśli fragment sygnału zostanie uznany jako poprawny
interwał VPC, będzie uwzględniony w obliczaniu TO i TS.
Warunki pozwalające na uznanie fragmentu jako poprawny VPC określono poniżej.
\item Dla potencjalnego interwału VPC zostaje wykonane obliczenie długości przedziału referencyjnego 
(średnia z  5 interwałów poprzedzających VPC).
\item Potencjalny RR z VPC musi być 
krótszy o przynajmniej 20 procent względem przedziału referencyjnego
\item Interwał występujący po potencjalnym interwale z VPC musi spełniać 
warunki pauzy kompensacyjnej (jego długość musi być większa o przynajmniej 10 procent od przedziału referencyjnego).
\item Pośród 5 interwałów przed VPC
oraz 20 interwałów po pauzie kompensacyjnej nie mogą znajdować się
interwały o długościach
mniejszych od 300 ms oraz większych od 2000 ms.
\item W ciągu 5 interwałów przed VPC
oraz 20 interwałów po pauzie kompensacyjnej musi występować
równy rytm  -- sąsiadujące interwały muszą różnić się o mniej
niż 200 ms oraz ich różnica względem przedziału referencyjnego musi być mniejsza
od 20 procent.
\item Jeśli interwał spełnia założone warunki, zostanie 
wykorzystany do obliczania parametrów TS oraz TO.
\item Na koniec sprawdzana jest liczba wykrytych 
PVC -- jeśli wykryto mniej niż 5, sygnał nie może być efektywnie 
wykorzystany, ponieważ nie zawiera odpowiedniej ilości
danych diagnostycznych. 
\end{itemize}
\textbf{Obliczanie TO}\\
\begin{figure}[h!]
\centering
%\includegraphics[width=15cm]{include/RRy.eps}
\caption{Schemat prezentujący umiejscowienie interwałów wykorzystywanych przy wyliczeniu parametru TS. 
\textit{ Źródło: http://www.h-r-t.org, Working Group of Biological Signal Analyses, Technische Universität München} }
\end{figure}
Dla poprawnie zidentyfikowanych VPC wyznaczany jest początek turbulencji.
Wykorzystywany jest przy tym wzór:
\begin{figure}[h!]
\centering
%\includegraphics[width=15cm]{include/ONSETequation.eps}
\caption{Wzór pozwalający na wyliczenie TO. }
\end{figure}
\noindent RR-2 i RR-1 to długości interwałów występujących bezpośrednio przed VPC,
a RR1 i RR2 to długości interwałów występujących bezpośrednio po pauzie 
kompensacyjnej.
Wstępnie TO określa się dla każdego pojedynczego VPC a następnie uśrednia się wszystkie wartości TO z analizowanego zapisu EKG.
\textbf{Obliczanie TS}\\
Algorytm obliczania turbulence slope:
\begin{itemize}
\item Na podstawie wyznaczonych odcinków zawierających poprawne VPC wyznaczany jest 
uśredniony tachogram.
\item Poddanych analizie zostaje 20 interwałów uśrednionego tachogramu występujących po pauzie kompensacyjnej. 
\item Spośród pierwszych dwudziestu interwałów po pauzie kompensacyjnej dla każdych kolejnych 5 interwałów wyznaczana
jest prosa regresji (brane są kolejno interwały 1-5, 2-6, 3-7 itd).
\item Spośród wyznaczonych prostych regresji zostaje wybrana ta ze współczynnikiem regresji 
o największej dodatniej wartości. 
\end{itemize}
Prostą regresji obliczamy wykorzystując współrzędne x, y pięciu kolejnych punktów należących
do uśrednionego tachogramu. Przy wykorzystaniu współrzędnych punktów tworzone są macierze
$A$ i $Y$.
\begin{equation}
A =  \left[
 \begin{array}{ccc}
   x_{1} & 1 \\
   x_{2} & 1 \\
   x_{3} & 1 \\ 
   x_{4} & 1 \\
   x_{5} & 1 \\ 
 \end{array}
\right]
\end{equation}

\begin{equation}
Y =  \left[
 \begin{array}{ccc}
   y_{1} \\
   y_{2} \\
   y_{3} \\ 
   y_{4} \\
   y_{5} \\ 
 \end{array}
\right]
\end{equation}
\\
Następnie, korzystając z równania  $(A' A)C = (A'Y)$ wyliczana jest macierz $C$.
Pierwszy element wynikowej macierzy $C$ to współczynnik kierunkowy prostej \textbf{a} (dla równania $y = ax + b$).
\subsubsection{Opis wykorzystanych w programie klas i funkcji}
\label{sec:hrt:results}

Poniżej opisano najważniejsze klasy, funkcje i zmienne wykorzystane w module:
   
Dołączona do projektu klasa \textbf{Matrix} jest klasą wykorzystywaną do wykonywania obliczeń. Dzięki niej można
definiować obiekty klasy Matrix reprezentujące dwuwymiarowe macierze oraz wykonywać na nich najważniejsze obliczenia.
Klasa Matrix została wykorzystana przy obliczaniu Turbulence Slope.

Klasa  \textbf{HRTAnalyzer} jest główną klasą modułu, dziedziczy po klasie HRTModule. Metody klasy HRTAnalyzer 
pozwalają na wyliczenie parametrów HRT.

Klasa \textbf{ECGHRT} jest klasą wykorzystywaną do przechowywania danych uzyskanych podczas pracy  -- przechowuje
najważniejsze obiekty będące wynikiem pracy HRTAnalyzera  -- parametry Turbulence Slope i Turbulence Onset, 
liczbę wykrytych VPC, oraz inne parametry wykorzystywane do rysowania.

Najważniejsze zmienne klasy ECGHRT:
\begin{itemize} 
\item Dane wykorzystywane przy rysowaniu wykresu: \textbf{QVector<QPointF> rr} wektor reprezentujący uśredniony tachogram, 
\textbf{QLine ts} -- linia reprezentująca Turbulence Slope.
\item Wyświetlane wartości:\textbf{int vpcCounter}  -- ilość znalezionych VPC, \textbf{double TO} --  średnia wartość Turbulence Onset,\textbf{double TS}  -- średnia wartość Turbulence Slope.
\end{itemize}


Funkcje i zmienne klasy HRTAnalyzer:
\begin{itemize}
\item \textbf{double* RRs}  -- zmienna wykorzystywana do tymczasowego przechowywania sygnału.
\item \textbf{void calculateHrtParams(double *signal, int size, ECGHRT \&)} -- funkcja, która zleca obliczenie parametrów HRT dla sygnału wejściowego i zapisuje je w wynikowym obiekcie. Parametry: signal to wskaźnik do tablicy interwałów RR, size to rozmiar tablicy, trzeci argument to referencja do obiektu przechowującego wyniki.
\item \textbf{vector<int> findVpcOnsets(double *signal, int size)} -- funkcja wyszukująca VPC w sygnale, zwracająca wektor 
z numerami próbek uznanymi za poprawne VPC, które mogą być wykorzystane przy działaniu programu. Parametry: signal to wskaźnik do tablicy interwałów RR, size to rozmiar tablicy.
\item \textbf{double* calculateAvgTach(double *signal, vector<int> vpc\_list) } -- funkcja obliczająca uśredniony tachogram.
Parametry: signal to wskaźnik do tablicy interwałów RR, vpc\_list to wektor indeksów oznaczajacych wystąpienie VPC.
\item \textbf{double calculateTO(double * signal, int size, vector<int> vpc\_list) }  -- funkcja obliczająca parametr TO,
Parametry: signal to wskaźnik do tablicy interwałów RR, size to rozmiar tablicy, vpc\_list to wektor indeksów oznaczających wystąpienie VPC.
\item \textbf{void calculateTS(double * signal, int size, vector<int> vpc\_list, double* avgTach, double to, ECGHRT\&)} -- funkcja obliczająca parametr TS. Parametry: signal to wskaźnik do tablicy interwałów RR, size to rozmiar tablicy, vpc\_list to wektor indeksów oznaczających wystąpienie VPC.
\end{itemize}

\subsubsection{Testy}
\label{sec:hrt:tests}

\begin{figure}[h!]
\centering
\includegraphics[width=15cm]{include/image.eps}
\caption{Zrzut ekranu prezentujący ostateczny wygląd modułu} 
\end{figure}

\begin{figure}[h!]
\centering
\includegraphics[width=15cm]{include/skrin1.eps}
\caption{Zrzut ekranu prezentujący ostateczny wygląd modułu} 
\end{figure}

\begin{figure}[h!]
\centering
\includegraphics[width=15cm]{include/skrin2.eps}
\caption{Zrzut ekranu prezentujący ostateczny wygląd modułu} 
\end{figure}

Dane do pierwszego etapu testów zostały uzyskane z bazy European ST-T Database (link: \url{http://www.physionet.org/}).
Początkowo testy wykonywano za pomocą prototypu aplikacji w Matlabie  -- testy działania aplikacji zakończyły się sukcesem 
(na wykresie wyświetlały się wartości TO i TS, wartości parametrów były zgodne z oczekiwaniami opartymi o analizę
wyglądu wykresów). 
Dla sygnału e0113 otrzymano wartości: turbulence slope -1.962, turbulence onset 6.03.
Dla sygnału e0123 otrzymano wartości: turbulence slope -0.957, turbulence onset 3.82.
Dla sygnału e0104 otrzymano wartości: turbulence slope -2.891, turbulence onset 8.272.
Niestety, w związku z faktem, że niektóre moduły, od których moduł HRT był zależny nie funkcjonowały
poprawnie na 13 godzin przed terminem oddania, poniechano dalszych bezowocnych prób dokonania bardziej miarodajnych
rodzajów testów na aplikacji. Nie mniej jednak, wcześniej dokonano podstawowej analizy działania modułu.
Testy polegały na ocenieniu poprawności rozwiązania w oparciu o analizę wyglądu wykresów. Dzięki wykresowi możliwa była
wzrokowa ocena poprawności otrzymywanych wyników (np. sprawdzenie, czy prosta regresji pokrywa się z najbardziej stromym zboczem
na wykresie, ocena wartości TO w odniesieniu do różnic długości). Testy nie wykazały uchybień. 


\bibliographystyle{plain}
\bibliography{report}
\end{document}
